# Re: Installation of Rb-GSL setup config problems

**URL:** <https://www.ruby-forum.com/t/re-installation-of-rb-gsl-setup-config-problems/94598>\
**Category:** Ruby\
**Created:** [May 8, 2006, 7:15pm UTC](https://www.ruby-forum.com/t/re-installation-of-rb-gsl-setup-config-problems/94598 "2006-05-08T19:15:58Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![unknown](https://www.ruby-forum.com/letter_avatar_proxy/v4/letter/u/4bbf92/32.png) [@unknown](https://www.ruby-forum.com/u/unknown)\
**Post date:** [May 8, 2006, 7:15pm UTC](https://www.ruby-forum.com/t/re-installation-of-rb-gsl-setup-config-problems/94598/1 "2006-05-08T19:15:58Z")

</div>

Dear Ara,

well thank you very much again.  
Actually, I am trying to find parameter values for statistical  
distributions  
from data, i.e.,  
you have, say, a normal distribution, and some measured data that are  
distributed  
according to it, but you don’t know the mean and variance parameters.  
Then, one can solve some equations called maximum likelihood equations  
for  
that.  
I was writing a Ruby script to do that. I wanted it to be able to find  
parameters for  
several well-known distributions. One of these is a lesser well-known  
one,  
the  
Cauchy distribution, which is sometimes used to model random events  
where  
extremely high values have a higher probability than if the data were  
distributed  
according to the normal distribution. Some people claim that, e.g.,  
natural  
disasters  
generate insurance claims that are distributed according to it or that  
the  
changes  
in stock market values are distributed like that ( when a crash  
occurs, and  
your money is lost, you know the famous (normal distribution)  
random-walk  
pricing models are wrong).  
This distribution has two parameters and I took a set of combined  
likelihood  
equations out of a book, which are solved numerically.

Now, thanks to you, I can solve them. Below is the code.

Best regards,

Axel

* * *

require “gsl”  
include Math  
include MultiMin

class Hash  
def ml\_fit  
fitted\_params={}

# fit distribution parameters from data using maximum likelihood or

method

# of moments

# self is a hash with keys ‘distribution’, and values any of the

following:

# 

‘normal’,‘bernoulli’,‘binomial’,‘uniform’,‘gamma’,‘exponential’,‘cauchy’

# ‘data’ - an Array with the data measured.

n=self[‘data’].length  
if self.has\_key?(‘distribution’)==false  
raise ‘No distribution specified’  
else  
if self[‘distribution’].downcase==‘normal’ or  
self[‘distribution’].downcase==‘gaussian’  
mu\_val=0.0  
self[‘data’].each{|x| mu\_val=mu\_val+x}  
mu\_val=mu\_val/self[‘data’].length  
fitted\_params.store(‘mu\_fitted’,mu\_val)  
sigma\_val=0.0  
self[‘data’].each{|x| sigma\_val=sigma\_val+(x-mu\_val)\*(x-mu\_val)}  
sigma\_val=sigma\_val/self[‘data’].length  
fitted\_params.store(‘sigma\_square\_fitted’,sigma\_val)  
end

if self[‘distribution’].downcase==‘bernoulli’

# p(x)=p^x\*(1-p)^(1-x), if x \in {0,1}, p(x)=0 otherwise

p\_fitted=0.0  
self[‘data’].each{|x| p\_fitted=p\_fitted+x}  
p\_fitted=p\_fitted/self[‘data’].length  
fitted\_params.store(‘p\_fitted’,p\_fitted)  
end  
if self[‘distribution’].downcase==‘binomial’

# needs parameter n !

# p(x)=\choose{n,x} p^x (1-p)^(n-x), if x \in {0,…,n}

if self.has\_key?(‘params’)==false  
raise(‘parameter n missing’)  
end  
p\_fitted=0.0  
self[‘data’].each{|x| p\_fitted=p\_fitted+x}  
p\_fitted=p\_fitted/self[‘data’].length/self[‘params’][‘n’]  
fitted\_params.store(‘p\_fitted’,p\_fitted)  
end  
if self[‘distribution’].downcase==‘poisson’

# p(x)=\frac{\lambada^x}{x!} exp(-\lambda), if x \in {0,1,…}

lambda\_fitted=0.0  
self[‘data’].each{|x| lambda\_fitted=lambda\_fitted+x}  
p\_fitted=p\_fitted/self[‘data’].length  
fitted\_params.store(‘lambda\_fitted’,lambda\_fitted)  
end  
if self[‘distribution’].downcase==‘uniform’

# p(x)=1/(b-a), if x \in [a,b]

a\_fitted=self[‘data’].max.to\_f  
b\_fitted=self[‘data’].min.to\_f  
fitted\_params.store(‘a\_fitted’,a\_fitted)  
fitted\_params.store(‘b\_fitted’,b\_fitted)  
end  
if self[‘distribution’].downcase==‘gamma’

# p\_{\alpha,\lambda}(x)=0, if x \leq 0;

# p\_{\alpha,\lambda}(x)=\frac{\lambda^{\alpha}}{\Gamma\_{\alpha}}

x^{\alpha-1}exp(-\lambda\*x), if x \>;

# use method of moments to estimate params

m\_1\_fitted=0.0  
self[‘data’].each{|x| m\_1\_fitted=m\_1\_fitted+x}  
m\_1\_fitted=m\_1\_fitted/self[‘data’].length

m\_2\_fitted=0.0  
self[‘data’].each{|x| m\_2\_fitted=m\_2\_fitted+x_x}  
m\_2\_fitted=m\_1\_fitted/self[‘data’].length  
alpha\_fitted=(m\_2\_fitted_m\_2\_fitted)/(m\_2\_fitted-m\_1\_fitted_m\_1\_fitted)  
lambda\_fitted=m\_1\_fitted/(m\_2\_fitted-m\_1\_fitted_m\_1\_fitted)  
fitted\_params.store(‘alpha\_fitted’,alpha\_fitted)  
fitted\_params.store(‘lambda\_fitted’,lambda\_fitted)  
end  
if self[‘distribution’].downcase==‘exponential’

# p(x)=\lambda_exp(-\lambda_x)

lambda\_fitted=0.0  
self[‘data’].each{|x| lambda\_fitted=lambda\_fitted+x}  
lambda\_fitted=lambda\_fitted/self[‘data’].length  
if lambda\_fitted==0  
raise(‘error fitting exponential distribution parameter: lambda is  
infinite’)  
else  
lambda\_fitted=1/lambda\_fitted  
fitted\_params.store(‘lambda\_fitted’,lambda\_fitted)  
end  
end  
if self[‘distribution’].downcase==‘cauchy’

# p(x,\alpha,\beta)=\frac{\beta}{\pi\*(\beta^2+(x-\alpha)^2))

# likelihood equations for parameters are

# 

# \frac{\partial L}{\partial \alpha}=0 \<=\> \sum\_{i=1}^{n}

\frac{x\_j-\alpha}{\beta^2+(x\_j-\alpha)^2}=0

# 

# \frac{\partial L}{\partial \beta}=0 \<=\> \sum\_{i=1}^{n}

\frac{\beta^2}{\beta^2+(x\_j-\alpha)^2}=n/2

# ( see Johnson, Kotz, and Balakrishnan, (1994), Continuous Univariate

Distributions,

# Volumes I and II, 2nd. Ed., John Wiley and Sons, chapter 16 for

details )

# Solve these equations using Nelder-Mead multidimensional

minimization

# do this several times to choose a good minimum from local minima

alpha\_v=[]  
beta\_v=[]

# now find the parameters, repeating the numerical minimization

repetitions=2000 # but stop after finding 5 parameter sets (in which  
the  
minimum search converged) to suitable values (beta\>0)  
temp\_min\_val=100000000.0  
best\_alpha=false  
best\_beta=false  
repetitions.times do  
fmin\_val,x,st=self.cauchy\_minimize  
if st==0 and x[1]\>0  
alpha\_v\<\<x[0]  
beta\_v\<\<x[1]  
end

if fmin\_val.abs\<temp\_min\_val.abs and x[1]\>0  
best\_alpha=x[0]  
best\_beta=x[1]  
temp\_min\_val=fmin\_val  
end

# stop if five good parameter values are found

if alpha\_v.length\>4  
break  
end  
if fmin\_val\<temp\_min\_val and x[1]\>0  
best\_alpha=x[0]  
best\_beta=x[1]  
end  
end  
=begin  
p ‘got back x’  
p alpha\_v  
p beta\_v  
p ‘best values’  
p best\_alpha  
p best\_beta  
=end  
fitted\_params.store(‘alpha\_fitted’,best\_alpha)  
fitted\_params.store(‘beta\_fitted’,best\_beta)  
end  
end  
return fitted\_params  
end

def cauchy\_minimize  
n=self[‘data’].length.to\_f  
np=self[‘data’].length  
values=self[‘data’]  
param\_fn=Proc.new{|dist\_params,data|

# params are measured values

val=0.0  
for x\_k in data  
val=val+(x\_k-dist\_params[0])  
/(dist\_params[1]_dist\_params[1]+(x\_k-dist\_params[0])_(x\_k-dist\_params[0]))  
end  
val=val.abs  
for x\_k in data  
val=val+(dist\_params[1]\*dist\_params[1])/(dist\_params[1]_dist\_params[1]+(x\_k-di  
st\_params[0])_(x\_k-dist\_params[0]))-n/2  
end  
val=val.abs  
}  
my\_func=Function.alloc(param\_fn,np)  
my\_func.set\_params(self[‘data’])  
x = Vector.alloc([rand]\*self[‘data’].length)  
ss = Vector.alloc(np)  
ss.set\_all(1.0)  
minimizer = FMinimizer.alloc(“nmsimplex”, np)  
minimizer.set(my\_func, x, ss)  
iter = 0  
begin  
iter += 1  
status = minimizer.iterate()  
status = minimizer.test\_size(1e-2)  
if status == GSL::SUCCESS

# puts(“converged to minimum at”)

end  
x = minimizer.x

# printf("%5d ", iter);

for i in 0…np do

# printf("%10.3e ", x[i])

end

# printf(“f() = %7.3f size = %.3f\n”, minimizer.fval, minimizer.size);

end while status == GSL::CONTINUE and iter \< 100  
return minimizer.fval,x,status  
end  
end

# test it

=begin  
a={‘distribution’,‘cauchy’,‘data’,[-1,1,2,-2]} # usually needs more  
data,  
of course, to get a good approximation  
params=a.ml\_fit  
p params  
=end