Fitting a Gaussian Curve

[Solabs]

dear all,

I built a simple model and tested it with a set of known values. I thought I could share the results here to get some feedback and hopefully get it done better and faster.

The model tries to fit a circular Gaussian, also known as von Mises function to a set of data points. Don't get confused, this is not a circular probability function, it is just a circular Gaussian curve with 3 parameters, the precision kappa, the difference between top and bottom, amp, and the offset of the curve i.e. the position on the y-axis where the minimum of the curve is located.

The data organized as clusters, so that each of the M clusters contain T such vM profiles. Thus each cluster has an hyperparameter, which dictates a prior on the profiles it contains. 

von Mises function is defined in the functions{} block... wrt performance, could one vectorize this operation, or wouldn't it matter any way?

between line 52-60, where I define the hyperparameters. should these lines be included in the for loop block?

I took the pain and generated some artificial data to compare the fitted parameters to their true values. I used 125 clusters that resulted by combining the following values... For each profile I added some random normal noise independently.

amp   =  [-5 -2.5 0 2.5 5];

offset  =  [-5 -2.5 0 2.5 5];

kappa =  [-4.6   -2.6  -0.7   1.3    3.2];

The figure depicts this relationship between fitted cluster parameters and ground truth. The second row shows the joint error matrix after marginalizing for the 3rd parameter.

I see that the offset and amplitude parameters are nicely fit and the values are close to the data  The kappa parameter is another story, there are lots of points all away from the identity line. But the joint error matrix shows that these are originating mostly from the case where the amplitude parameter is 0, which makes the kappa value undetermined. Otherwise, I am a bit surprised about the sigma_noise parameter which is always undershot by the fitting procedure even though it has a rather non-informative prior (line 56).

cheers,

s.

functions {

        #Function to fit to dependent variables. This is a vonMises function scaled to span [0 1],

        #centered on 0 degree.

        #amp is the difference between peak and thorough

        #offset is the height of the minimum value.

        #kappa is the precision parameter.

        #no free parameter for the location of the peak, it is always at 0 degrees.

        real vonMises(real x, real amp, real kappa, real offset) {

                #how to vectorize this?, cos does only take real values not real[].

                real y;

                y <- amp*((exp(kappa*cos(x))-exp(-kappa))/(exp(kappa)-exp(-kappa)))+offset;

                return y;

        }

}

data {

        #the data consists of T columns of X data points, which are to be modelled using the vM function

        #Columns originates M different groups and within each group there T profiles.

        #Each participant is present in each group:

        #Nth profile in each group comes from the same participant, but this is ignored here. 

        int<lower=0>                     X;              #number of independent variables.

        vector[X]                        x;              #independent variables in degrees.

        int<lower=0>                     T;              #total number of profiles.

        real                             y[X,T];         #dependent variable, each column is one vMises profile.

        int<lower=0>                     M;              #number of groups.

        int                              m[T];           #indicator variable for origin of y (specifies the recording modality).

}

parameters {

        #group-level parameters:

        real<lower=0>                    sigma_y[M];     #noise variance for each group

        real                             mu_amp[M];      #average amplitude for the group

        real                             mu_offset[M];   #average offset for the group

        real<lower=0,upper=400>          mu_kappa[M];    #average precision for the group, values bigger than 400 leads to NaN output from vonMises function above.

        real<lower=0>                    sigma_amp[M];   #variance for amplitude for group

        real<lower=0>                    sigma_offset[M];#variance for offset for the group

        real<lower=0>                    sigma_kappa[M]; #variance for kappa

        #profile level parameters

        real                             amp[T];         #amplitude for single profiles 

        real                             offset[T];      #offset    for single profiles

        real<lower=0,upper=400>          kappa[T];       #precision for single profiles

}

transformed parameters {

}

model {

                #group level (does this need to be inside the for loop?)        

                sigma_kappa   ~ cauchy(0,1);#what alternative are there for these parameters?

                sigma_amp     ~ cauchy(0,1);

                sigma_offset  ~ cauchy(0,1);

                sigma_y       ~ cauchy(0,10);

                mu_kappa      ~ lognormal(0,1.5);

                mu_amp        ~ cauchy(0,10);

                mu_offset     ~ cauchy(0,10);

                for (i in 1:T) {#run over profiles      

                        #parameters of single profile are distributed according to group-level parameters.

                        #could this be vectorized?

                        offset[i]           ~ normal(mu_offset[m[i]]         , sigma_offset[m[i]]);

                        kappa[i]            ~ lognormal( log(mu_kappa[m[i]]) , sigma_kappa[m[i]]);

                        amp[i]              ~ normal(mu_amp[m[i]]            , sigma_amp[m[i]]);

                        for (xi in 1:X){#run over the dependent variable                

                                y[xi,i] ~ normal(  vonMises(x[xi],amp[i],kappa[i],offset[i]) , sigma_y[m[i]] );

        }}}

[Mike Lawrence]

Stan already has von_mises_lpdf(), so you don't really need to have your own function. Indeed, I believe von_mises_lpdf() accepts vectors for its arguments. To implement your 3-parameter model, you'd just be doing "amp*von_mises_lpdf(y|offset,kappa)"

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Mike Lawrence]

Oh, and you'll want to restrict the range of the offset parameter, else you'll get multiple modes as -pi is the same as +pi

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[Solabs]

hi,

hmm I am not trying to fit a probability function. also in your description the amp parameter would be hard to interpret as it would have different scaling effect depending on kappa. in my vM function amp sets the peak value.

best,

[Mike Lawrence]

Doesn't matter. von_mises_lpdf just gives the height of the pdf at the point y given the parameters offset & kappa. Actually, it gives the log of the height, so you'll want to do "amp*exp(von_mises_lpdf(y|offset,kappa)"

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[Solabs]

offset is actually the global height of the curve on the y-axis, not its position on the x-axis (which is fixed at 0 degrees)...

On Monday, December 5, 2016 at 11:14:26 AM UTC+1, Solabs wrote:

[Solabs]

it does matter. 

it will certainly fit something, but the parameter amp will depend on kappa. just see this way: if kappa is 0.1 or 1, same amplitude value will give different peak values.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Mike Lawrence]

Sorry, yes, I read your initial post too quickly and mis-guessed your parameterization. You should be able to do:

    exp(von_mises_lpdf( y | 0,kappa))/exp(von_mises_lpdf( 0 | 0,kappa))*(amp-offset)+offset

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

In general, we can't map a circular statistic to R
in a way that preserves continuity. I've been meaning to
try using our unit_vector parameterization for points on
a circle, which overaparameterizes and constrains, but
does have locality of the right kind.

You can't learn much from point estimates in Bayesian analysis;
you want to measure coverage of intervals to see if the parameters
are being recovered properly. In general, there can be biases
if you build a hierarchical model and the data aren't generated
from such a model.

- Bob

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Solabs]

thx Bob for the conceptual remarks, made me thing about my analysis.