@Carl Love 

Thanks.

And can you explain the code you gave me before?

Thanks.

@Carl Love 

f:= E-> exp(-E/kT)/A;
eq1:=f(0)=2.143;
eq2:=f(1)=1.484;
eq3:=f(2)=0.989;
solve({eq1,eq2,eq3},{A,kT});

But I can't see the root.

@Carl Love 

So, how can I get the answer?

@Preben Alsholm 

I have another  problem.

The code is:

f:=E->B/(A*exp(E/kT));
eq1:=f(0)=2.143;
eq2:=f(1)=1.484;
eq3:=f(2)=0.989;
solve({eq1,eq2,eq3},{A,kT,B});
plot(eval(f(E),%),E=0..9);

But it  runs error.

@Preben Alsholm 

Thank you very much.

@Preben Alsholm 

I have problem.

res:=Statistics:-NonlinearFit(f(E),M,E);

Can you explain the code ?

Thank you very much.

@Preben Alsholm 

Thank you very much.

@Carl Love 

Thank you for your answer. But I can't understand your code.

Can you explain your the code approximately ?

 Can you help me ? Thanks

@Carl Love 

Sorry I didn't notice it. I corrected it

@Carl Love 

Thank you.  But I want to use the same random sequences  to compile  three different statistics.

I am sorry I forget to say.

ex:

[1, 6, 0, 0, 2, 0], [1, 1, 1, 4, 1, 1], [0, 2, 0, 0, 7, 0], [4, 3, 1, 0, 1, 0], [2, 0, 0, 0, 0, 7]

Maxwell Boltzmann statistics

[0 = 13/5, 1 = 8/5, 2 = 3/5, 3 = 1/5, 4 = 2/5, 6 = 1/5, 7 = 2/5]

Bose-Einstein statistics

The sequence will be [1, 6, 0, 0, 2, 0], [1, 1, 1, 4, 1, 1], [0, 2, 0, 0, 7, 0], [4, 3, 1, 0, 1, 0]

or [1, 6, 0, 0, 2, 0], [1, 1, 1, 4, 1, 1], [4, 3, 1, 0, 1, 0], [2, 0, 0, 0, 0, 7].

[0 = 9/4, 1 = 8/4, 2 = 2/4, 3 = 1/4, 4 = 2/4, 6 = 1/4, 7 = 1/4]

fermi-dirac statistics

The sequence will be [4, 3, 1, 0, 1, 0] .

[0 =2/1, 1=2/1 ,3=1/1 , 4=1/1 ]

@Carl Love 

Thank you very much for your help. But that you give the code is inclined to the theory.

My subjet is random.  I want to use the code that you provided before to vertify the three types of the theory about statistics.

RandomCompositions:= proc(n::posint, k::posint)
local
C,
Compositions:= [seq(C-~1, C= combinat:-composition(n+k, k))],
Rand:= rand(1..nops(Compositions))
;
()-> Compositions[Rand()]
end proc:

R:= RandomCompositions(9,6):
n:= 100000:
S:= 'R()' $ n:
T:= map(lhs=rhs/n, Statistics:-Tally(op~([S]))):
T:= [{T[]}[]]:

You gave me the code can vertify the Maxwell-Boltzmann distribution.

I want to use the samplers up to the point " ten thousand or more" random sequences to vertify  Bose-Einstein condensation and Fermi-Dirac statistics. Thanks.

@Carl Love 

You are right. The website is wrong.

The Maxwell-Boltzmann distribution sampling times are 2002.

But you just list five sequences. I don't know why.

I am very confuced about the sampling. How can you set sampling times are 2002?

You said the computation is exact---not based on random sampling. But latter you said that does either the exact computation or provides a random sampler depending on an input option.

What's the difference from random sampling and random sampler

@Carl Love 

Oh,I forgot it.I am so sorry.Thank you very much.

But I still have many problem.

How can I set  the sampling frequency?

Population:= select(p-> nops(p)=k, combinat:-partition(n+k, n+1));

I think that is Bose-Einstein statistics.

What's different from you gave me before?

The Bose-Einstein energy 4 is 0.269.

But the  result is 0.3076. It has  much inaccuracy. 

This is the Average number. 

http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/disfdx.html#c2

Thank you very much.

 But I can't plot. 

I am sorry I confused the Bose-Einstein condensation and Fermi-Dirac statistics.