Skip to main content
Mathematica

Return to Question

added 452 characters in body
Source Link
YoAn
YoAn
  • 35
  • 6

EDIT: Following @JimB comment, below is the code I used to try to estimate parameters with fse using data that were directly generated from it:

data = RandomVariate[fse[1, 0.1, 2, 0.01], 100];
FindDistributionParameters[data, fse[so, l, n, sd] /. sd -> 0.01, 
  ParameterEstimator -> {"MaximumLikelihood", Method -> "NMaximize", MaxIterations -> 500, AccuracyGoal -> 5}];

It also runs indefinitely

EDIT: Following @JimB comment, below is the code I used to try to estimate parameters with fse using data that were directly generated from it:

data = RandomVariate[fse[1, 0.1, 2, 0.01], 100];
FindDistributionParameters[data, fse[so, l, n, sd] /. sd -> 0.01, 
  ParameterEstimator -> {"MaximumLikelihood", Method -> "NMaximize", MaxIterations -> 500, AccuracyGoal -> 5}];

It also runs indefinitely

Correct title spelling
Link
mikado
mikado
  • 17.1k
  • 2
  • 21
  • 55

How to include measurement error in FindDsitributionParametersFindDistributionParameters (using maximum likelihood)?

Source Link
YoAn
YoAn
  • 35
  • 6

How to include measurement error in FindDsitributionParameters (using maximum likelihood)?

I am trying to fit a non-standard PDF to data using FindDistributionParameters. My custom distribution is defined like this:

fs = Function[{so,l, n},
   TransformedDistribution[
    so - l/2 c,
    c \[Distributed] 
     NoncentralChiSquareDistribution[n, 2 so/l],
     Assumptions -> so > 0 && n > 0 && l > 0
    ]
   ];

My data corresponds to measurements of a variable on N individuals. For each individual I have multiple measurements to take into account the measurement error. My data is similar to this :

SeedRandom[1];
trueValue = RandomVariate[fs[so, l, n] /. {so -> 1, l -> 0.1, n -> 2}, 100];(* true value for each individual, which I don't know*)
errorSD = RandomVariate[ChiSquareDistribution[3], 100]/10; (*standard deviations of the measurment error for each individual*)

fullMeasurement = RandomVariate[NormalDistribution[#[[1]], #[[2]]], 5] & /@ Transpose[{trueValue, errorSD}];
myData = {Mean[#], StandardDeviation[#]} & /@ fullMeasurement(*Simulated data set similar to my real data*)

As a first approach, I used FindDistributionParameters on the mean of the measurements.

p = FindDistributionParameters[myData[[All,1]],
   fs[so, l, n],
   ParameterEstimator -> {"MaximumLikelihood", Method -> "NMaximize",
     MaxIterations -> 500, AccuracyGoal -> 5}];

Show[
 Histogram[myData[[All,1]], 15, "ProbabilityDensity"],
 Plot[PDF[fs[so,l, n] /. p, x], {x, -2, 1}, PlotRange -> All, PlotStyle -> Thick]
 ]

Yellow : histogram of simulated data. Blue line : PDF of the custom distribution with the estimated parameters

However with this approach I am not using the information of the standard deviation of the measurements. As a second approach, assuming the error is normally distributed, I tried to make a convolution of my custom distribution and a normal distribution. I first tried to get the PDF of the convolution using either Convolve or directly integrating with Integrate but it did not work.

Then tried to use TransformedDistribution to create the convolution and use FindDistributionParameters:

fse = Function[{so, l, n, sd},
   TransformedDistribution[
    so - l/2 c + e,
    {c \[Distributed] NoncentralChiSquareDistribution[n, 2 so/l],
     e \[Distributed] NormalDistribution[0, sd]},
    Assumptions -> so > 0 && n > 0 && l > 0 && sd > 0
    ]
   ];

FindDistributionParameters[data,
  fse[so, l, n, sd] /. sd -> Mean[data[[All,2]]],
  ParameterEstimator -> {"MaximumLikelihood", Method -> "NMaximize", 
    MaxIterations -> 500, AccuracyGoal -> 5}];

However it runs indefinitely, and this approach is using only the mean of the standard deviations and not the value for each individual. Is their a better (and more efficient) way to include the measurement error in the estimation method ?