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

stat: implement Pareto #309

Merged
merged 1 commit into from Dec 8, 2017

Conversation

Projects
None yet
3 participants
@kczimm
Copy link
Contributor

commented Nov 21, 2017

An implementation of the Pareto (Type I) distribution without the Rand random sampling and the Quantile methods.

@kczimm kczimm force-pushed the kczimm:stat-pareto branch 2 times, most recently from c6c394e to 2f986f1 Nov 21, 2017

@kczimm

This comment has been minimized.

Copy link
Contributor Author

commented Nov 21, 2017

Added the Rand method and associated tests. Still do not have the Quantile method implemented though.

// with support above the scale parameter.
//
// The density function is given by
// \frac{\alpha x_m^{\alpha}}{x^{\alpha+1}} for x >= x_m

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member

We use something that is almost but not quite TeX for these formulae. In this case I'd not use \frac but rather a / and parens. Also, indent the formula by one space and finish with a dot.

// The density function is given by
// \frac{\alpha x_m^{\alpha}}{x^{\alpha+1}} for x >= x_m
//
// For more information, see https://en.wikipedia.org/wiki/Pareto_distribution

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member

Final full stop.

This comment has been minimized.

Copy link
@kczimm

kczimm Nov 21, 2017

Author Contributor

I don't know the proper lingo yet. What does final full stop mean? I put a trailing period here, but I'm not sure if that's what you meant.

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member

I don't speak ASE, "full stop" == "period".

// Xm is the scale parameter. Must be greater than 0.
Xm float64

// Alpha is the shape parameter. Must be greater than 0.

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member

"Alpha must be greater than 0."

//
// For more information, see https://en.wikipedia.org/wiki/Pareto_distribution
type Pareto struct {
// Xm is the scale parameter. Must be greater than 0.

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member

"Xm must be greater than 0."

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member

Sorry, here and below, only change the second sentence. The first sentence was fine and needs to be there.

if p.Alpha <= 4 {
return 0
}
return 6 * (math.Pow(p.Alpha, 3) + math.Pow(p.Alpha, 2) - 6*p.Alpha - 2) / p.Alpha / (p.Alpha - 3) / (p.Alpha - 4)

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member

s/math.Pow(p.Alpha, 3)/p.Alphap.Alphap.Alpha/
s/math.Pow(p.Alpha, 2)/p.Alpha*p.Alpha/
s|p.Alpha / (p.Alpha - 3) / (p.Alpha - 4)|(p.Alpha * (p.Alpha - 3) * (p.Alpha - 4))|


// Mean returns the mean of the probability distribution.
func (p Pareto) Mean() float64 {
if p.Alpha > 1 {

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member
if p.Alpha <= 1 {
    return math.Inf(1)
}
return p.Alpha * p.Xm / (p.Alpha - 1)

// Prob computes the value of the probability density function at x.
func (p Pareto) Prob(x float64) float64 {
if x < p.Xm {

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member

Usually we return math.Exp(p.LogProb(x)).

This comment has been minimized.

Copy link
@kczimm

kczimm Nov 21, 2017

Author Contributor

done

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member

Omit the condition, you already check for x < p.Xm in LogProb and return -Inf in that case; math.Exp(math.Inf(-1)) returns 0.

This comment has been minimized.

Copy link
@kczimm

kczimm Nov 21, 2017

Author Contributor

Ah, so obvious. Thank you. Done.

if p.Alpha <= 2 {
return math.Inf(1)
}
return p.Xm * p.Xm * p.Alpha / math.Pow(p.Alpha-1, 2) / (p.Alpha - 2)

This comment has been minimized.

Copy link
@kortschak

kortschak Nov 21, 2017

Member
am1 := p.Alpha-1
return p.Xm * p.Xm * p.Alpha / (am1 * am1 * (p.Alpha - 2))

@kczimm kczimm force-pushed the kczimm:stat-pareto branch from 628180b to 51bb509 Nov 21, 2017

// with support above the scale parameter.
//
// The density function is given by
// (\alpha x_m^{\alpha})/(x^{\alpha+1}) for x >= x_m.

This comment has been minimized.

Copy link
@btracey

btracey Dec 6, 2017

Member

Could you swap \alpha with α

if x < p.Xm {
return 0
}
return 1 - math.Pow(p.Xm/x, p.Alpha)

This comment has been minimized.

Copy link
@btracey

btracey Dec 6, 2017

Member

This should return 1 - p.Survival(x)

AUTHORS Outdated
@@ -38,6 +38,7 @@ Joseph Watson <jtwatson@linux-consulting.us>
Josh Wilson <josh.craig.wilson@gmail.com>
Julien Roland <juroland@gmail.com>
Kent English <kent.english@gmail.com>
Kevin C. Zimmerman <kevinczimmerman@gmail.com>

This comment has been minimized.

Copy link
@btracey

btracey Dec 6, 2017

Member

We may not need these if the other one gets in first.

}
return p.Alpha * p.Xm / (p.Alpha - 1)
}

This comment has been minimized.

Copy link
@btracey

btracey Dec 6, 2017

Member

Median? It's easy here.

This comment has been minimized.

Copy link
@kczimm

kczimm Dec 7, 2017

Author Contributor

done.


// Rand returns a random sample drawn from the distribution.
func (p Pareto) Rand() float64 {
return p.Xm / math.Pow(1-Uniform{Min: 0, Max: 1}.Rand(), 1/p.Alpha)

This comment has been minimized.

Copy link
@btracey

btracey Dec 6, 2017

Member
  1. This doesn't need to construct a distribution, that's built into the *rand.Rand
  2. It does need to use the rand.Rand to generate the variable if it's non-nil
  3. This should work in log space instead
    math.Exp(math.Log(p.Xm) - p.Alpha * rnd.ExpFloat64()) (I think that's the right sign)
if x < p.Xm {
return 1
}
return math.Pow(p.Xm/x, p.Alpha)

This comment has been minimized.

Copy link
@btracey

btracey Dec 6, 2017

Member

I would also do this in log space and convert it back. math.Exp(p.Alpha*( math.Log(p.Xm) - math.Log(x)))

}
return 1 - math.Pow(p.Xm/x, p.Alpha)
}

This comment has been minimized.

Copy link
@btracey

btracey Dec 6, 2017

Member

Entropy?

for _, test := range []struct {
x, xm, alpha, want float64
}{
{0, 1, 1, 0},

This comment has been minimized.

Copy link
@btracey

btracey Dec 6, 2017

Member

Try some non-integer values in here too please.

This comment has been minimized.

Copy link
@kczimm

kczimm Dec 7, 2017

Author Contributor

done.

@kczimm kczimm force-pushed the kczimm:stat-pareto branch from 09282a6 to 51bb509 Dec 7, 2017

@btracey
Copy link
Member

left a comment

LGTM with the entropy change.


// Entropy returns the differential entropy of the distribution.
func (p Pareto) Entropy() float64 {
return math.Log((p.Xm / p.Alpha) * math.Exp(1+1/p.Alpha))

This comment has been minimized.

Copy link
@btracey

btracey Dec 7, 2017

Member

math.Log(p.Xm) - math.Log(p.Alpha) + (1 + 1/p.Alpha)

@kczimm kczimm force-pushed the kczimm:stat-pareto branch from 8a958fc to 6562f38 Dec 7, 2017

@kczimm

This comment has been minimized.

Copy link
Contributor Author

commented Dec 7, 2017

Fixed entropy and squashed.

@btracey

btracey approved these changes Dec 7, 2017

@btracey btracey merged commit 153d483 into gonum:master Dec 8, 2017

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.