mmcdara

6214 Reputation

17 Badges

7 years, 313 days

MaplePrimes Activity


These are replies submitted by mmcdara

@Axel Vogt 

Thnk you Axel, this really helps.

@jalal 

Let's say we have some sample S, for instance S := Statistics:-Sample(Normal(0, 1), 10^3).
What is it that does not suit you in the command Statistics:-Histogram(S, opts), where opts is any sequence of options described in the help page?

Here are a few examples of histograms of samples drawn from a continuous or discrete random variable.
The last example is the histogram of the table you provided in your original question.

restart:

with(Statistics):


Histogram of a sample from a continuous random variable.

S := Sample(Uniform(-5, 5), 10^3):

h0 := Histogram(S):
h1 := Histogram(S, style=polygon):
h2 := Histogram(S, maxbins=10):
h3 := Histogram(S, frequencyscale=absolute):
h4 := Histogram(S, binbounds=[seq(-5+k, k=0..10)]):
h5 := Histogram(S, frequencyscale=absolute, transparency=0.8, color=red):

plots:-display(Matrix(2, 3, [seq(h||k, k=0..5)])):


Histogram of a sample from a discrete random variable.

S := Sample(DiscreteUniform(-5, 5), 10^3):


h0 := Histogram(S):
h1 := Histogram(S, discrete=true):
h2 := Histogram(
        S
        , discrete=true
        , axis[1]=[tickmarks=[$(-5..5)]]
        , thickness=20
        , style=polygon
        , view=[-5.5..5.5, default]
      ):

plots:-display(Matrix(1, 3, [h0, h1, h2])):


Histogram of  your sample.

data :=[ [[50, 60], 20], [[60, 70], 25], [[70, 80], 15], [[80, 90], 30], [[90, 100], 10] ]

[[[50, 60], 20], [[60, 70], 25], [[70, 80], 15], [[80, 90], 30], [[90, 100], 10]]

(1)

S := [seq((add(data[k][1])/2)$data[k][2], k=1..numelems(data))]:

Histogram(
  S
  , frequencyscale=absolute
  , discrete=true
  , axis[1]=[tickmarks=[seq(k=cat("[", k-5, ", ", k+5, "["), k in [seq(55..95, 10)])]]
  , thickness=40
  , style=polygon
  , view=[50..100, default]
);

 

 

Download histograms.mw

Please let me/us know if you have a particular need. For instance provide data (or the way they are constructed) and an image of the histogram you would obtain.

First thing: what you represent IS NOT an Histogram (look the corresponding help pages).

Second point "[You are] looking to construct a histogram in an aesthetic manner with bins. The Histogram command does not provide an optimal output. Any suggestions?".... what does aesthetic  mean to you? Are you aware the what is aesthetic for you is a subjective judgement and may seem awfullto someone else?

So, begin to provide an example of what you feel is an aesthetic histogram for you and we will be capable to give you a solif answer.

Nevertheless the attached file contains a few examples of customized "Classes-Effectifs" tables.
Note that as I use Maple 2015 I don't know what your own version proposes to display a table like th one in your mw file.
So here hre a few examples of what can be done TallyInto.mw.

Last remark: If I'm not mistaken the lines

InsertContent(Worksheet(Group(Input( TAB )))):

must be written

InsertContent(Worksheet( TAB )):

since Maple 2018 (read help page DocumentTools:-Layout:-Cell)

Other usefil info can be found here Statistics:-TallyInto and Statistics:-FrequencyTable

@Axel Vogt 

Your two suggestions 

M := convert(Int(Dirac(u-v), [u, v]), piecewise, u);
W := convert(Int(Dirac(u+v-1), [u, v]), piecewise, u);

are very good.

The two following steps 

M-piecewise(min(u, v), u) ;
W-piecewise(max(u+v-1, 0), u) ;

leave me a little frustrated for they mean that we suspect that M could be equal to min(u,v) and that we verify this hypothesis.
I would have prefer to apply successive transformations to finaly get (Oh miracle!) the expression min(u, v).

I nevertheless vote up for the the frst idea.

Remove the characters "&", "(", ")" in the name of your file: it can be loaded (page not found)

@fabs 

Here is the edited file which includes the case X=x1*tauderivatives.mw
In fact the trick works whatever the function of x1 and tau.

@Carl Love 

Thanks Carl, I understand better now.

@Carl Love 

Than you Carl.
This is almost what I did just after havng read @vv's answer.

restart
phi := x -> (x^(-theta)-1)/theta:
ihp := u -> solve(phi(x)=u, x):
invfunc[phi] := ihp:

C := (u, v) -> (phi@@(-1))(phi(u)+phi(v)):
simplify(C(u, v)) assuming theta >= -1, theta <> 0;
                                         /    1  \
                                         |- -----|
                                         \  theta/
             /      (-theta)    (-theta)\         
             \-1 + u         + v        /         

The only slight difference (but maybe of great importance???) is that I did not use unprotect(invfunc).
I only checked out that this table indeed contained phi=ihp.

 

@vv 

Thank you vv.
I'm feeling sorry for being so stupid because I had looked at the invfunc help page but obviously not carefully enough.

@acer 

Thank you for this suggestion and the related comments

I only knew about the operator-form calling sequence for plotting commands ot in the case of Optimization, where I use it on a regular basis.
I didn't know it could also be used int too.

And yes "the approach can break due to some overlooked extra evaluation which strips off a pair of such quotes": I unfortunately made this bad experience a lot of times.

@Preben Alsholm 

Thanks you Preben.
I had just edited my previous reply meanwhile.

@Preben Alsholm 

Sometimes I write my procedures the way you suggest, in particular when these procedures are used in a plot command.
But I also noted that these definitions both gave the same correct result:

MyProc1 := proc(u)
  if not u::realcons then return 'procname'(_passed) end if; 
  fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = u, x)
end proc:

plot(MyProc1(u), u=0.1..0.9):


MyProc2 := proc(u)   
  fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = u, x) 
end proc: 

plot('MyProc2'(u), u=0.1..0.9):


Son lazzy as I am, I often use the second definition.
This is why I used it in my question.

So, it seems that defining the procedure as in MyProc1 is "safer' than in MyProc2. Am I right?
Do you systematically recommend proceeding as in MyProc1 ?


By the way: the computational time comes from the fact that 

fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = q, x)

is equal to -infinity if q=0 and +infinity if q=1. In practice I would integreate within the range 0.01..0.99 for instance.

For the record here is the method I used to use to numerically evaluate J2 before I thought of using evalf/Int instead (last block in the attached file).
I suppose that a judicious choice of the method in evalf/Int can provide the result fastly then the default method used in J3 ?

 

restart:

J2 := proc()
  local z:
   z := proc(q1, q2)  
     local x1, x2;
     if not [q1,q2]::list(realcons) then return 'procname'(_passed) end if;
     x1 := fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = q1, x):
     x2 := fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = q2, x):
     0.2652582384*exp( 2*x1*x2 - x1^2 - x2^2 )
  end proc:
  evalf(Int(z(q1, q2), q1=0.01..0.99, q2=0.01..0.99))
end proc:

CodeTools:-Usage( J2() );

Warning,  computation interrupted

 

J3 := proc()
  local z:
   z := proc(q1, q2)  
     local x1, x2;
     uses Statistics:
     if not [q1,q2]::list(realcons) then return 'procname'(_passed) end if;
     x1 := Quantile(Normal(0, 1), q1, numeric):
     x2 := Quantile(Normal(0, 1), q2, numeric):
     exp( 2*x1*x2 - x1^2 - x2^2 )
  end proc:
  evalf(Int(z(q1, q2), q1=0.01..0.99, q2=0.01..0.99))
end proc:

CodeTools:-Usage( J3() );

memory used=0.76GiB, alloc change=44.01MiB, cpu time=14.31s, real time=14.55s, gc time=689.80ms

 

.4449972351

(1)


A more efficient way

StringTools:-FormatTime("%H:%M:%S");

h := 0.001:
P := [seq(0.01..0.99, h)]:
Q := CodeTools:-Usage( Statistics:-Quantile~(Normal(0, 1), P, numeric) ):

S:= 0:
for q1 in Q do
  for q2 in Q do
    S := S + exp( 2*q1*q2 - q1^2 - q2^2 )
  end do:
end do:
S * h^2;

StringTools:-FormatTime("%H:%M:%S");
 

"14:39:25"

 

memory used=16.68MiB, alloc change=0 bytes, cpu time=291.00ms, real time=292.00ms, gc time=0ns

 

HFloat(0.44516891241912115)

 

"14:39:29"

(2)

 


 

Download Computational_times.mw

@Ronan 

Writing 

{clr::string:="Blue"}

in the argument list of a procedure means that clr is to be considered as an optional argument whose type is string and default value "Blue".
If you do not mention clr when you execute the procedure spread2, clr will take inside this procedure the default value "Blue".
If you want the procedure uses "Red" instead of "Blue", you have to set ecplivitely clr="Red" when you call it.

Optional arguments are intensively used in Maple: you may see them as an annoyance for they impose writing clr="Red" instead of simply "Red", but their advantage is that the order ofthe optional parameters doesn't matter (so you can write indifferently clr="Red", prnt=false, or prnt=false, clr="Red").
I greatly appreciate this feature, but it remains a personal opinion.

@acer 

Thank you acer.
Your

eval(y, g=''f'')

suits me perfectly.

By the way I sent at @Preben Alsholm a file where the real problem is described, in case you would be interested.

@Preben Alsholm 

Thanks four answer.

"why not simply start with f=F?" simply because this example was notional (maybe oversimplified).
Here is a more detailed example which is closer to the reality (note that writing it I realized there was a simple way [end of the file] to answer my own question!)

 

restart

with(LinearAlgebra):
with(Statistics):

phi := unapply(CDF(Normal(0, 1), x), x);

proc (x) options operator, arrow; 1/2+(1/2)*erf((1/2)*x*2^(1/2)) end proc

(1)

Sigma := Matrix(2$2, [1, rho, rho, 1]):
X := Vector(2, symbol=x):

phi__Sigma := exp(-1/2*(X^+ . Sigma^(-1) . X)) / (2*Pi*sqrt(Determinant(Sigma)));

(1/2)*exp(-(1/2)*(-x[1]/(rho^2-1)+x[2]*rho/(rho^2-1))*x[1]-(1/2)*(x[1]*rho/(rho^2-1)-x[2]/(rho^2-1))*x[2])/(Pi*(-rho^2+1)^(1/2))

(2)

# Let psi the CDF of some continuous random variable V.
# We define C y

C__psi := eval(phi__Sigma, [seq(x[i]=((phi@@(-1))@cdf)(q[i]), i=1..2)]);
 

proc (x) options operator, arrow; psi(x) end proc

 

(1/2)*exp(-(1/2)*(-(phi@@(-1))(psi(q[1]))/(rho^2-1)+(phi@@(-1))(psi(q[2]))*rho/(rho^2-1))*(phi@@(-1))(psi(q[1]))-(1/2)*((phi@@(-1))(psi(q[1]))*rho/(rho^2-1)-(phi@@(-1))(psi(q[2]))/(rho^2-1))*(phi@@(-1))(psi(q[2])))/(Pi*(-rho^2+1)^(1/2))

(3)

# In the specific case V is a standard gaussian random variable we have cdf(V) = phi

C__psi := eval(phi__Sigma, [seq(x[i]=((phi@@(-1))@phi)(q[i]), i=1..2)]);

 

(1/2)*exp(-(1/2)*(-q[1]/(rho^2-1)+q[2]*rho/(rho^2-1))*q[1]-(1/2)*(q[1]*rho/(rho^2-1)-q[2]/(rho^2-1))*q[2])/(Pi*(-rho^2+1)^(1/2))

(4)

# But

cdf := x -> phi(x):
C__psi := eval(phi__Sigma, [seq(x[i]=((phi@@(-1))@cdf)(q[i]), i=1..2)]);

(1/2)*exp(-(1/2)*(-(phi@@(-1))(1/2+(1/2)*erf((1/2)*q[1]*2^(1/2)))/(rho^2-1)+(phi@@(-1))(1/2+(1/2)*erf((1/2)*q[2]*2^(1/2)))*rho/(rho^2-1))*(phi@@(-1))(1/2+(1/2)*erf((1/2)*q[1]*2^(1/2)))-(1/2)*((phi@@(-1))(1/2+(1/2)*erf((1/2)*q[1]*2^(1/2)))*rho/(rho^2-1)-(phi@@(-1))(1/2+(1/2)*erf((1/2)*q[2]*2^(1/2)))/(rho^2-1))*(phi@@(-1))(1/2+(1/2)*erf((1/2)*q[2]*2^(1/2))))/(Pi*(-rho^2+1)^(1/2))

(5)

 

Writing this more realistic example revealed me how to proceed in a simpler way

 

C:= psi ->  eval(phi__Sigma, [seq(x[i]=((phi@@(-1))@psi)(q[i]), i=1..2)]):

C(psi);
C(phi);

(1/2)*exp(-(1/2)*(-(phi@@(-1))(psi(q[1]))/(rho^2-1)+(phi@@(-1))(psi(q[2]))*rho/(rho^2-1))*(phi@@(-1))(psi(q[1]))-(1/2)*((phi@@(-1))(psi(q[1]))*rho/(rho^2-1)-(phi@@(-1))(psi(q[2]))/(rho^2-1))*(phi@@(-1))(psi(q[2])))/(Pi*(-rho^2+1)^(1/2))

 

(1/2)*exp(-(1/2)*(-q[1]/(rho^2-1)+q[2]*rho/(rho^2-1))*q[1]-(1/2)*(q[1]*rho/(rho^2-1)-q[2]/(rho^2-1))*q[2])/(Pi*(-rho^2+1)^(1/2))

(6)

 


 

Download RealProblem.mw

 

2 3 4 5 6 7 8 Last Page 4 of 126