Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Nonuniform distributions for Morris sampling returns inf #515

Open
willu47 opened this issue Jul 12, 2022 · 11 comments
Open

Nonuniform distributions for Morris sampling returns inf #515

willu47 opened this issue Jul 12, 2022 · 11 comments

Comments

@willu47
Copy link
Member

willu47 commented Jul 12, 2022

Thanks to @jyangfsu for the issue submission.

When the parameters follows norm, lognorm and truncnorm distributions, nonuniform_scale_samples function would returns inf.
As suggested by Saltelli et al. (2010), this can be avoid by cutting the tails of , for example, the normal distributions, at quantiles 5 and 95%.
Code in current nonuniform_scale_samples function:

elif dists[i] == 'norm':
    if b2 <= 0:
        raise ValueError("""Normal distribution: stdev must be > 0""")
    else:
        conv_params[:, i] = sp.stats.norm.ppf(params[:, i], loc=b1, scale=b2))

Suggest modifying to:

elif dists[i] == 'norm':
    if b2 <= 0:
        raise ValueError('''Normal distribution: stdev must be > 0''')
    else:
        conv_params[:, i] = scipy.stats.truncnorm.ppf(
            params[:, i], scipy.stats.norm.ppf(0.05, loc=b1, scale=b2), scipy.stats.norm.ppf(0.95, loc=b1, scale=b2), loc=b1, scale=b2)
@ConnectedSystems
Copy link
Member

Hi @jyangfsu

Is the reference below the one you are referring to? I had a quick scan over it but could not see exactly where your suggestion is mentioned so I think you are referring to a different paper.

Saltelli, A., P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola (2010).
"Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index."
Computer Physics Communications, 181(2):259-270,
doi:10.1016/j.cpc.2009.09.018.

@jyangfsu
Copy link

jyangfsu commented Oct 11, 2022 via email

@sitadrost
Copy link

Sorry, bit late to the party, but I just installed SALib (version 1.4.5, which is the newest available through conda) because I was looking for a package to do Morris screening in Python. As my variables are normally distributed, I stumbled into this issue straight off.
The suggested modification is to use a truncated normal distribution, which sounds sensible, but I don't think the suggested code is quite correct.
From the truncnorm docs: The standard form of this distribution is a standard normal truncated to the range [a, b], where a and b are user-provided shape parameters. The parameter loc shifts the mean of the underlying normal distribution, and scale controls the standard deviation of the underlying normal, but a and b are still defined with respect to the standard normal.
So I think the modified code should simply be:

elif dists[i] == 'norm':
    if b2 <= 0:
        raise ValueError('''Normal distribution: stdev must be > 0''')
    else:
        conv_params[:, i] = scipy.stats.truncnorm.ppf(
            params[:, i], 0.05, 0.95, loc=b1, scale=b2)

Or, to capture 95% of the normal distribution in the truncated version (instead of 90%, like in the code above):

elif dists[i] == 'norm':
    if b2 <= 0:
        raise ValueError('''Normal distribution: stdev must be > 0''')
    else:
        conv_params[:, i] = scipy.stats.truncnorm.ppf(
            params[:, i], 0.025, 0.975, loc=b1, scale=b2)

@tupui
Copy link
Member

tupui commented Oct 11, 2022

(version 1.4.5, which is the newest available through conda)

Thanks for pointing this out. It's surprising to me, the system should have noticed the update on PyPi and make a PR to update to 1.4.6. I will make a PR to update.

@ConnectedSystems
Copy link
Member

Thanks @tupui - I made a mental note to check on the conda feedstock after the 1.4.6 release but obviously it had left my mind when I woke in the morning.

@tupui
Copy link
Member

tupui commented Oct 11, 2022

@sitadrost
Copy link

Woops, I just noticed that the modified code that I proposed last week doesn't give the results you'd expect: it returns values between the mean (b1) and the mean + 1x standard deviation (b2), instead of values ranging from (b1 - 2 x b2) to (b1 + 2 x b2), like you'd expect for a normal distribution. Not quite sure yet where exactly this goes wrong, have you looked at this in the meantime?

@ConnectedSystems
Copy link
Member

Thanks @sitadrost , I'll try to make some time to have a look this weekend.

@sitadrost
Copy link

Ah, of course, so silly of me: truncnorm works with a standard normal distribution, so the code to capture 95% of the desired normal distribution (mean plus/minus 2 x standard deviation) should be:

elif dists[i] == 'norm':
    if b2 <= 0:
        raise ValueError('''Normal distribution: stdev must be > 0''')
    else:
        conv_params[:, i] = scipy.stats.truncnorm.ppf(
            params[:, i], -2, 2, loc=b1, scale=b2)

@ConnectedSystems
Copy link
Member

ConnectedSystems commented Dec 3, 2022

Apologies for the long silence on this issue, short on time these days.

I suspect that the suggested code will have implications for other sampling approaches, as the scaling is not used for just the Morris method. The suggested changes cause the test for the Oakley function to fail, as an example. This test uses LHS and expects all parameters to be normally distributed.

I think in this case, would it be advisable to just direct users to specify truncnorm rather than norm with whatever tail cut-off is desired.

That said, I'm not sure this is the "best" solution. For example, where to signpost not to use norm for normally distributed factors in the case of Morris? I could put in a warning/error for users of the ProblemSpec interface, and assume those using the functions directly are power users who know what they're doing...

Any thoughts here @willu47 @tupui ?

@khairulislam
Copy link

Hi, I am facing the same issue in the latest version. Has it been solved?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

6 participants