39350_en-US

liliconcom

Peak of a probability distribution

Hi

I have a dataset which fits best to the Johnson Su distribution and I was wondering if there is a way to estimate the value of the peak of this distribution. It has just one peak, which is pretty different from the median or the mean. Calling this "mode" might be the best possible solution here but how do I estimate this point?

Re: Peak of a probability distribution

Can the profile option be extended for data distribution that is grouped? For N groups, need to get the mode corresponding to distribution of each of the items in the group. Or should I subset by group everytime and then use the profiler option. I used the col groupby and then fitted but looks like it is giving a single value for entire dataset and not for each group member.

Re: Peak of a probability distribution

You have a solution using the profiler that will work. As for the scripts not working for you, this article highlights some of the changes made to the Maximize and Minimize functions. Thanks to the JMP Tech Support team to pointing it out to me!

https://community.jmp.com/t5/JMPer-Cable/Minimize-and-Maximize-Functions-in-JSL/ba-p/36355

 

Re: Peak of a probability distribution

Laura's method is neater! And far less error prone. It's always best to let JMP do the work for you if you can figure out a way.

Re: Peak of a probability distribution

One of the things that I love about JMP is that it often gives you many ways to solve a problem.  All of the previous suggestions for finding the mode of the fitted Johnson SU density for the data are great, but here is one more option that uses the Profiler.

 

Once you have used Distribution to fit the Johnson Su distribution to your data, you can choose Save Density Formula from the Fitted Johnson Su menu.  This saves the Johnson Su Density formula to a column in your data table.

Next you can open the Profiler by choosing Graph>>Profiler.  Cast your density formula as the Y, Prediction Formula and click OK.  Next turn on the Desirability Functions from Prediction Profiler menu.  After that you choose Maximize Desirability from the Prediction Profiler menu.  This will find the mode (maximum) of your Johnson Su density for this data.

 

 

 

 

 

Re: Peak of a probability distribution

Thank you so much!

Re: Peak of a probability distribution

I am not having any problem with the data in this column. Here is a new version that avoids having to manually enter the parameter estimates, borrowing from Ian's solution:

Names Default to Here( 1 );

dist = Current Data Table() << Distribution(
 Stack( 1 ),
 Continuous Distribution(
  Column( :T10_100 ),
  Horizontal Layout( 1 ),
  Vertical( 0 ),
  Fit Distribution( Johnson Su )
 )
);

estimates = Report( dist )["Fitted Johnson Su"][NumberColBox(1)] << Get As Matrix;

dist expr = Substitute(
 Expr(
  Johnson Su Density(
   q,  // quantile
   ggg, // gamma or first shape parameter
   ddd, // delta or second shape parameter
   ttt, // theta or location parameter
   sss  // sigma or scale parameter
  )
 ),
 Expr( ggg ), estimates[1],
 Expr( ddd ), estimates[2],
 Expr( ttt ), estimates[3],
 Expr( sss ), estimates[4]
);

peak height = Maximize( dist expr, { q( 0.01, 70 ) } );

difference = Minimize( peak height - dist expr, { q( 0.01, 70 ) } );

Print( "The peak location is " || Char( q ) );

I get this result with your data:

 

So it is not the data or the script, it must be JMP (installation).

Re: Peak of a probability distribution

The optimisation doesn't work in the sample data table also. Also, I used Ian's script which is completely independent of my data. So there must be something wrong with the installation then, as you pointed it out. Should I reinstall JMP now?

I'm attaching the data column anyway.

Re: Peak of a probability distribution

I don't see anything wrong with Ian's script or mine so that leaves the data table or the JMP installation. Let's start with the data table. Can you share it? If so, please identify the data column that you used. Actually, we only need the data column that is causing the problem. You can delete the rest if that is better for you.

Thanks!

Re: Peak of a probability distribution

It didn't work for your script and for my data set. Same with Ian's script, it says

"Optimization failed: Failed: Cannot Decrease Objective Function".

So I just went ahead with bill's solution, but I consider it as only a temporary fix for my problem.

Re: Peak of a probability distribution

Are you saying that my solution script didn't work with the original example that you posted (as it did for me) or that it didn't work with a new data set? If it is the latter case, can you share the data?

Did you try Ian's script?

Re: Peak of a probability distribution

Mark's method is neat. Here's an alternative that uses the function directly:

NamesDefaultToHere(1);
ClearLog();

// Make some representative data
dt = NewTable("Johnson Su", NewColumn("Data", Numeric, Continuous, Formula(Random Johnson Su( 5, 2, 1, 1 ))));
dt << addRows(50);

// Fit the data
dist = dt << Distribution(Continuous Distribution( Column( :Data ), Fit Distribution( Johnson Su ) ));
// Get the fitted parameters
params = Report(dist)[NumberColBox(3)] << get;
gamma = params[1]; // Shape
delta = params[2]; // Shape
theta = params[3]; // Location
sigma = params[4]; // Scale

// http://www.mathwave.com/articles/johnson_su_distribution.html
//   f(x) = delta/(lambda * Sqrt(2*pi()) * Sqrt(z^2 + 1)) * exp( -1/2*(gamma + delta * ln(z + Sqrt(z^2+1)))^2 );
//   where:
//   z = (x - eta)/lambda;
//  and:
//   delta > 0 lambda > 0
jsu = Expr(delta/(lambda * Sqrt(2*pi()) * Sqrt(z^2 + 1)) * exp( -1/2*(gamma + delta * ln(z + Sqrt(z^2+1)))^2 ) );
SubstituteInto(jsu, Expr(z), Expr((x - eta)/lambda));
// Use the fitted parameter values
// Note the change in parameter names! 'lambda' -> 'sigma' and 'eta' -> 'theta'
SubstituteInto(jsu,
 Expr(gamma), Eval(gamma),
 Expr(delta), Eval(delta),
 Expr(eta), Eval(theta),
 Expr(lambda), Eval(sigma)
);
// Get upper and lower bounds for the optimisation (either side of the (assumed) single peak)
vals = Column(dt, "Data") << getValues;
low = Quantile(0.25, vals);
high = Quantile(0.75, vals);
// Do the maximization
Maximize(jsu, {x(low, high)}, << Tolerance(10^-10), << showDetails(true));
Print(x);

('By eye' it seems to work, but no testing of course).

Re: Peak of a probability distribution

Thank you for the response. But somehow it doesn't work for me. It says

"Optimization failed: Failed: Maximum Iterations Exceeded"

I'm using JMP 12, maybe the optimizer has a limited functionality.

Re: Peak of a probability distribution

Just to add on.

Mark's way is really the best way to get the peak value, but for a quick way to get an estimated value you can fit the curve then use the JMP Crosshairs tool found in the Tools menu or on your menu bar.  Click the tool then put the crosshairs symbol on the spot you feel best represents the peak and a value will pop up.  If you are using JMP 13 you will get a magnifier with the Crosshair tool that will help you get close. 

HTH

Re: Peak of a probability distribution

Unfortunately, the Distribution platform does not provide the mode for an analytical function. I searched but could not find the formula for the mode of a Johnson distribution either.

Another approach is to use the JMP numerical optimizer to search for the solution. My approach is to use the Maximize() function to find the maximum density of the fitted distribution, which corresponds to the peak, and then use the Minimize() function to find the quantile the produces a probability closest to the maximum. Here is the solution in steps to match what you are trying to do.

1. I opened the Injection Molding example from the Sample Data Directory.
2. I opened the Distribution platform on the Shrinkage data column.
3. I fit the Johnson SU distribution to obtain the Johnson SU parameter estimates.
4. I wrote a script with two steps outlined above to find the peak location.

Here is the fitted distribution:

.The script is:

peak height = Maximize(
 Johnson Su Density(
  q,   // quantile
  -25.83224, // gamma or first shape parameter
  1.5352821, // delta or second shape parameter
  -3.912839, // theta or location parameter
  2.5166e-6 // sigma or scale parameter
 ),
 { q( 0.01, 70 ) }
);

difference = Minimize(
 peak height - Johnson Su Density(
  q,   // quantile
  -25.83224, // gamma or first shape parameter
  1.5352821, // delta or second shape parameter
  -3.912839, // theta or location parameter
  2.5166e-6 // sigma or scale parameter
 ),
 { q( 0.01, 70 ) }
);

Print( "The peak location is " || Char( q ) );

Notice that the estimates are copied into the function call, which is used twice.

 

Here is the result in the Log:

Don't be alarmed by all of the messages - it worked just fine! You can see the result at the bottom.