curso9.mw

Capítulo 9  Probabilidade, Estatística e Análise de Dados 

Ref. Probabilidade e Estatística, M.R. Spiegel, J. Schiller, A. Srinivasan, terceira edição. 

Introdução 

O principal pacote para estatística é Statistics . Vamos também carregar outros pacotes que serão úteis em determinados pontos. 

> restart

> with(Statistics); -1

> with(LinearAlgebra); -1

> with(ListTools); -1

> with(plots); -1

> interface(rtablesize = 100); -1

A lista de comandos do pacote Statistics pode ser encontrada através do comando ? Statistics[Commands] . 

 

Problema. Um novo cassino começou a operar e proprôs os seguinte jogo. A banca tem três cartas: um ás e dois valetes. As três cartas são colocadas viradas para baixo de forma aleatória. Você deve escolher uma delas. Se a carta escolhida for um ás, você ganha US$10.000,00, se for um valete, você paga US$10.000,00. Porém, antes de abrir a carta escolhida, a banca abre uma das outras duas cartas e mostra onde está um valete, depois lhe pergunta, você quer continuar com a mesma carta originalmente escolhida ou prefere trocar para a outra carta que ainda está fechada. Qual opção abaixo é a mais correta. 

a) Você troca de carta pois a chance de ganhar um ás passa a ser de dois em três. 

b) Você não troca de carta pois é irrelevante para aumentar a chance de ganhar, que é um em três. 

c) Você troca de carta pois a chance de ganhar um ás passa a ser maior do que um em três, porém menor do que dois em três. 

Solução 

Este problema tem uma história muito interessante. Ele apareçeu pela primeira vez em uma carta de um estatístico e foi usado em um programa de TV pelo apresentador Monty Hall. O problema do Monty Hall ficou famoso quando foi usado por Marilyn vos Savant em uma coluna de perguntas em uma revista americana chamada Parade. A solução dada por Marilyn dizia que se você trocar de carta, a probabilidade de ganhar é 2/3 e se continuar com a mesma carta, a probabilidade de ganhar é 1/3. Em função desta resposta, Marilyn recebeu em torno de 10.000 cartas de leitores indignados, dizendo que ela estava errada e que não havia a menor necessidade de trocar de carta, pois isto não aumentava a chance de ganhar. Cerca de 1.000 cartas eram de matemáticos e em especial, o famoso matemático Paul Erdős também estava no grupo dos que achavam que não adiantava trocar de carta. Acontece que a Marilyn já tinha aparecido no Guiness como o recorde de melhor QI de todos os tempos. Depois destas informações, você prefere a opção a) ou b)? Parece que a opção c) é carta fora do baralho. Se você não trocar de carta, a probabilidade de ganhar é 1/3 independentemente do que o apresentador fizer. Agora vamos calcular a probabilidade de PERDER caso você opte pela estratégia de trocar de carta. Se você tiver escolhido o ás inicialmente  e trocar de carta, você perde. Se você tiver escolhido um valete, e trocar de carta você ganha, pois o outro valete foi revelado pela banca. Portanto, a probabilidade de PERDER neste caso é 1/3; a probabilidade de ganhar é 2/3. 

Variável Aleatória e Distribuição de Probabilidades 

Distribuição de probabilidades contínua 

Uma variável aleatória X é uma variável cujo valor é determinado de modo aleatório dentro de um espaço amostral. Cada valor de X tem uma probabilidade associada, que é descrita por uma distribuição de probabilidades. Uma distribuição de probabilidades pode ser discreta ou contínua dependendo da natureza do espaço amostral. A distribuição normal é uma distribuição de probabilidades contínua (o espaço amostral é o conjunto dos números reais) enquanto que a distribuição binomial é discreta (o espaço amostral é um subconjunto do números inteiros). A distribuição normal com média mu = 10 e desvio padrão sigma = 1 é dada por 

> randomize(); -1

> dn := PDF(Normal(10, 1), t)
`+`(`/`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(`+`(t, `-`(10)), 2))))))))), `*`(`^`(Pi, `/`(1, 2))))) (2.1.1)
 

cujo gráfico é 

> plot(dn, t = 6 .. 14)
Plot_2d
 

O comando PDF é um acrônimo de Probability Density Funtion e retorna a expressão analítica de uma distribuição de probabilidades. O comando Normal(mu, sigma) retorna a distribuição normal de média mu e desvio padrão sigma. A figura acima mostra o gráfico da distribuição normal de média 10 e desvio padrão 1. Uma distribuição de probabilidades deve sempre satisfazer as seguintes propriedades: (1) os valores da distribuição tem que estar no intervalo [0, 1] e (2) a área embaixo da curva deve ser 1. Note que 

> int(dn, t = `+`(`-`(infinity)) .. infinity)
1 (2.1.2)
 

Uma maneira alternativa de fazer o gráfico é usando diretamente a partir distribuição de probabilidades (sem construir a expressão analítica)  é usando o comando DensityPlot do pacote Statistics : 

> DensityPlot(Normal(10, 1), range = 6 .. 14, thickness = 2, color = blue)
Plot_2d
 

Uma variável aleatória é definida através do comando RandomVariable(dist) , cujo único argumento dist deve ser uma distribuição de probabilidades. Por exemplo, vamos definir uma variável aleatória com uma distribuição normal de média 10 e desvio padão 1: 

Para obter valores desta variável, devemos usar o comando Sample(X,n) onde X é a variável aleatória e n é o tamanho do vetor resposta. Por exemplo: 

> X := RandomVariable(Normal(10, 1)); -1

> V := Sample(X, 5)
V := rtable(1 .. 5, [9.95836578220813, 9.36584564756441, 10.2614261233570, 10.2653773533797, 9.84440640579403], subtype = Vector[row]) (2.1.3)
 

A média dos valores do vetor linha V (média amostral) deve ser aproximadamente 10, de fato 

> media := `/`(`*`(add(V[i], i = 1 .. Dimension(V))), `*`(Dimension(V)))
HFloat(9.93908426246066) (2.1.4)
 

No entanto, o valor correto da média associada à distribuição de probabilidades é calculado com o comando Mean : 

> Mean(X)
10 (2.1.5)
 

e o desvio padrão através do comando StandardDeviation : 

> StandardDeviation(X)
1 (2.1.6)
 

Note que o comando Sample retorna um vetor linha como resposta. O Maple não mostra as entradas se o vetor for grande, pois exite um limite de corte. Para aumentar o limite, é necessário usar o comando interface (rtablesize=1000), por exemplo, para mostrar até 1000 entradas. 

 

A lista das distribuições contínuas implementas no Maple é 

Beta Cauchy ChiSquare Erlang Error Exponential FRatio Gamma Gumbel InverseGaussian Laplace Logistic LogNormal Mawell Moyal NonCentralBeta NonCentralChiSquare NonCentralFRatio NonCentralStudentT Normal Pareto Power Rayleigh StudentT Triangular Uniform VonMises Weibull  

 

 

Exercício 

Suponha que uma distribuição de probabilidades contínua tem a forma f(x) = `/`(`*`(c), `*`(`+`(`*`(`^`(x, 4)), 1))) onde c é uma constante e `and`(`<`(`+`(`-`(infinity)), x), `<`(x, infinity)).  (a) Determine o valor de c. (b) Calcule a probabilidade de `*`(`^`(x, 2)) estar entre 1/3 e 1. 

 

Se f(x) é uma distribuição de probabilidades contínua, a função de distribuição acumulada é definida como 

F(x) = int(f(t), t = `+`(`-`(infinity)) .. x) 

Por exemplo, a função de distribuição acumulada da distribuição de probabilidades Normal(10, 1) é 

>

> fda := CDF(Normal(10, 1), x)
`+`(`/`(1, 2), `*`(`/`(1, 2), `*`(erf(`+`(`*`(`/`(1, 2), `*`(`+`(x, `-`(10)), `*`(`^`(2, `/`(1, 2)))))))))) (2.1.7)
 

cujo gráfico é 

> plot(fda, x = -1 .. 20)
Plot_2d
 

A função de distribuição acumulada sempre começa do 0 ( e termina no 1 (. Várias outras terminologias são usadas como função de distribuição cumulativa, função de distribuição, função densidade, entre outras. Em Inglês, usa-se o nome cumulative distribution function (CDF). 

Distribuição de probabilidades discreta 

Vamos ver um exemplo de uma variável aleatória com uma distribuição de probabilidades discreta. Vamos definir Y como uma variável aleatória com a distribuição binomial com n = 40 tentativas e probabilidade de sucesso p = `/`(1, 2): 

>

A distribuição binomial Binomial(n, p) tem com espaço amostral os números inteiros de 0 a n. Ela tem a seguinte interpretação. Suponha que temos uma moeda cuja probabilidade de dar CARA é p. Vamos supor que CARA é sucesso.  O valor Binomial(x, p) é a probabilidade de termos x CARAS em n tentativas. Por exemplo 

> Y := RandomVariable(Binomial(40, `/`(1, 2))); -1

> interface(rtablesize = 100)
100 (2.2.1)
 

> L := Sample(Y, 20)
L := rtable(1 .. 20, [20., 23., 18., 19., 24., 22., 22., 22., 26., 19., 24., 18., 23., 23., 17., 23., 26., 26., 22., 27.], subtype = Vector[row]) (2.2.2)
 

A lista L é interpretada da seguinte forma. Cada elemento da lista L é o número de CARAS em 40 jogadas da moeda. Note que a maioria dos números da lista L está perto de 20, pois a probabilidade de dar 20 CARAS em 40 jogadas é maior do que qualquer outro número de CARAS. A distribuição de probabilidades é encontrada usando o comando ProbabilityFunction(Y,x) , onde Y é a variável aleatória e x é um número inteiro entre 0 e n. O gráfico da distribuição é encontrado da seguinte maneira: 

> g1 := plot([seq([x, evalf(ProbabilityFunction(Y, x))], x = 0 .. 40)], style = point, symbol = diagonalcross, symbolsize = 20); -1; g1
g1 := plot([seq([x, evalf(ProbabilityFunction(Y, x))], x = 0 .. 40)], style = point, symbol = diagonalcross, symbolsize = 20); -1; g1
Plot_2d
 

A média e o desvio padrão de Y são 

> mu := Mean(Y)
20 (2.2.3)
 

> sigma := StandardDeviation(Y)
`*`(`^`(10, `/`(1, 2))) (2.2.4)
 

A distribuição normal e a distribuição binomial têm uma forte correlação. Vamos fazer o gráfico sobreposto destas distribuições: 

> g2 := plot(PDF(Normal(mu, sigma), i), i = 0 .. 40); -1

> display([g1, g2])
Plot_2d
 

Podemos ver que os pontos da distribuição binomial estão muito próximos da curva da distribuição normal, porém não exatamente. Por exemplo, para i = 20: 

> evalf(`+`(PDF(Normal(mu, sigma), 20), `-`(ProbabilityFunction(Y, 20))))
0.7859384e-3 (2.2.5)
 

A distribuição de probabilidades pode ser aproximada pelo histograma ( histogram ) de grandes amostras, por exemplo 

> Histogram(Sample(Y, 1000))
Plot_2d
 

A lista de distribuições discretas implementadas no Maple é 

Bernoulli Binomial DiscreteUniform EmpiricalDistribution Geometric Hypergeometric NegativeBinomial Poisson ProbabilityTable  

Exercícios 

1. (a) Gere uma amostra com 500 valores da variável Y (definida acima) e salve em uma lista L. (b) Calcule o número de ocorrências do número 20 na amostra L do ítem (a). (c) Seja N(i) o número de ocorrências do número i na lista L. Faça o gráfico dos valores de `+`(`*`(`/`(1, 500), `*`(N(i)))) para i de 0 a 40. (d) Qual é a relação do gráfico do ítem (c) e o gráfico da distribuição de probabilidades? Faça o histograma de L e descreva a relação entre o histograma o gráfico do ítem (c). 

 

2. Considere uma fábrica de garrafas. Suponha que uma máquina desta fábrica produz 5% de garrafas defeituosas. Suponha que selecionamos aleatoriamente uma amostra de 20 garrafas produzidas por esta máquina. (a) Qual é a probabilidade da amostra conter 3 garrafas defeituosas? (b) Qual é a probabilidade da amostra ter mais do que 3 garrafas defeituosas? 

 

3. Um laboratório está testando duas vacinas. A probabilidade da primeira vacina gerar reação alérgica é de 0.5% e da segunda é 2.8%. A primeira é testada em um grupo de 230 pessoas e a segunda em um grupo de 100 pessoas. Qual é a probabilidade de que, no total, 4 pessoas tenham reação alérgica? 

 

4. Suponha as notas de um exame nacional tem distribuição normal com média 7.0 e desvio padrão 1.2. Uma amostra aleatória de 100 estudantes é selecionada. Qual é a probabilidade de que a nota mais baixa do grupo dos 20% melhores seja menor do que 7.8? 

 

 

Se f(x) é uma distribuição de probabilidades discreta, a função distribuição acumulada (discreta) é definida como 

 

Por exemplo, a função distribuição acumulada da variável aleatória Y é 

>

> Y := RandomVariable(Binomial(40, `/`(1, 2))); -1

> F := unapply(sum(ProbabilityFunction(Y, t), t = 0 .. x), x)
(2.2.6)
 

O gráfico de F é 

> plot([seq([x, F(x)], x = 0 .. 40)], style = point)
Plot_2d
 

 

Gerando novas distribuições de probabilidades contínuas 

Para construir uma nova distribuição de probabilidades contínua (diferente das distribuições implementadas no pacote Statistics, veja Distributions ), devemos primeiro definir uma função no Maple que descreve a distribuição de probabilidades, por exemplo, seja a função f 

>

> f := proc (z) options operator, arrow; piecewise(`and`(`<=`(0, z), `<=`(z, 1)), 1, 0) end proc
proc (z) options operator, arrow; piecewise(`and`(`<=`(0, z), `<=`(z, 1)), 1, 0) end proc (2.3.1)
 

cujo gráfico é 

> plot(f(z), z = -1 .. 2)
Plot_2d
 

A função  f satisfaz as condições para ser uma distribuição de probabilidades, pois  f(z) está no intervalo [0,1] para todo z e a área embaixo da curva é 1, como podemos verificar com o seguinte comando 

> int(f(z), z = 0 .. infinity)
1 (2.3.2)
 

Para criar uma nova distribuição de probabilidades contínua, devemos usar o comando Distribution com a seguinte sintaxe 

> mydist := Distribution(PDF = f)
module () export Conditions, PDF; option Distribution, Continuous; end module (2.3.3)
 

A distribuição é um módulo ( module ). Qualquer outro resultado que não seja "module( ) option Distribuion, ..." indica que houve um erro. Agora podemos usar a distribuição mydist como se fosse qualquer outra distribuição de probabilidades contínua do Maple. Por exempo, podemos definir uma nova variável aleatória coma distribuição mydist e fazer uma amostragem: 

> Z := RandomVariable(mydist); -1

> Sample(Z, 5)
rtable(1 .. 5, [.833845868631608, .595316888630301, .834371984639605, 0.409843378372628e-1, .415837747258859], subtype = Vector[row]) (2.3.4)
 

Podemos calcular a média e o desvio padrão: 

> Mean(Z)
`/`(1, 2) (2.3.5)
 

> StandardDeviation(Z)
`+`(`*`(`/`(1, 6), `*`(`^`(3, `/`(1, 2))))) (2.3.6)
 

O desvio padrão pode ser confirmado através da definição 

> sqrt(int(`*`(`^`(`+`(z, `-`(`/`(1, 2))), 2), `*`(f(z))), z = `+`(`-`(infinity)) .. infinity))
`+`(`*`(`/`(1, 6), `*`(`^`(3, `/`(1, 2))))) (2.3.7)
 

 

Exercício 

Seja g(x) a função `*`(`*`(`+`(`*`(3, `*`(x))), `/`(1, 2)), `*`(`+`(1, `-`(`*`(`/`(1, 2), `*`(x)))))) no intervalo [0,2] e 0 fora deste intervalo. (a) Defina uma distribuição de probabilidades contínua com g(x). (b) Defina uma variável aleatória X associada e faça uma amostragem de tamanho 11. (c) Calcule a média e desvio padrão de X. (d) Seja Y uma variável aleatória com a distribuição de probabilidades g(x). Defina a variável aleatória Z = `+`(X, Y). Faça o gráfico da curva de distribuição de probabilidades de Z e mostre que a área embaixo da curva é igual a 1. 

Gerando novas distribuições de probabilidades discretas 

Para definir uma nova distribuição de probabilidades discreta, é necessário fornecer o espaço amostral e as probabilidades associadas a cada ponto do espaço amostral. Por exemplo, suponha que o espaço amostral seja infinity. Suponha que a probabilidade associada ao ponto n seja `/`(1, `*`(`^`(2, n))).  Temos que especificar quais são as probabilidades associadas a cada ponto através de uma função da seguinte maneira. 

>

> prob := proc (n) options operator, arrow; `/`(1, `*`(`^`(2, n))) end proc
proc (n) options operator, arrow; `/`(1, `*`(`^`(2, n))) end proc (2.4.1)
 

> seq(prob(n), n = 1 .. 10)
`/`(1, 2), `/`(1, 4), `/`(1, 8), `/`(1, 16), `/`(1, 32), `/`(1, 64), `/`(1, 128), `/`(1, 256), `/`(1, 512), `/`(1, 1024) (2.4.2)
 

Note que o valor de cada probabilidade está no intervalo [0,1] e além disso a soma das probabilidades é 

> sum(prob(n), n = 1 .. infinity)
1 (2.4.3)
 

De fato, toda distribuição de probabilidades discreta tem que satisfazer estas duas condições: (1) o valor de cada probabilidade deve estar no intervalo [0,1] e (2) a soma das probabilidades deve dar 1. Para definir a distribuição de probabilidade, devemos usar o comando Distribution com a seguinte sintaxe: 

onde Support = 1 .. infinity indica a os valores da variável n usada em DiscreteValueMap. No caso acima, a função que descreve os pontos do espaço amostral é proc (n) options operator, arrow; n end proc. Neste caso, n vai assumir os valores 1, 2, ..., infinity.  Vamos agora definir uma variável com esta nova distribuição de probabilidade e fazer uma amostragem 

> dist := Distribution(Type = discrete, Support = 1 .. infinity, DiscreteValueMap = (proc (n) options operator, arrow; n end proc), ProbabilityFunction = prob); -1

> Z := RandomVariable(dist); -1

> Sample(Z, 10)
rtable(1 .. 10, [2., 1., 3., 1., 2., 1., 1., 1., 1., 5.], subtype = Vector[row]) (2.4.4)
 

A média e o desvio padrão são 

> Mean(Z)
2 (2.4.5)
 

> StandardDeviation(Z)
`*`(`^`(2, `/`(1, 2))) (2.4.6)
 

 

Exercícios 

1. Suponha que o espaço amostral seja o conjunto dos pontos 1, 2, 3, ..., infinity. Seja f  uma distribuição de probabilidades discreta tal que a probabilidade associada ao ponto n do espaço amostral seja f(n) = `/`(`*`(c), `*`(`^`(3, n))). (a) Determine o valor de (b) Faça o gráfico da distribuição de probabilidades. (c) Faça o gráfico da função de probabilidades acumulada. (d) Calule P(`<=`(X, 6)) onde X é uma variavel aleatória com distribuição f. (e) Calule P(`>`(X, 3)). 

 

2. Suponha que o espaço amostral seja o conjunto de pontos `/`(1, 2), `/`(1, 4), `/`(1, 8), `/`(1, 16), `/`(1, 32), `/`(1, 64), () .. () , infinity. Defina uma distribuição de probabilidades discreta através do comando Distribution tal que a probabilidade associada ao ponto n do espaço amostral seja `/`(1, `*`(`^`(2, n))). Faça uma amostragem com 11 pontos, calcule a média e o desvio padrão de uma variável aleatória com esta distribuição de probabilidades e faça o gráfico da distribuição de probabilidades. 

Solução 

>

> pontos := proc (n) options operator, arrow; `^`(`/`(1, 2), n) end proc
proc (n) options operator, arrow; `^`(`/`(1, 2), n) end proc (2.4.1.1.1)
 

> dist := Distribution(Type = discrete, Support = 1 .. infinity, DiscreteValueMap = pontos, ProbabilityFunction = (proc (x) options operator, arrow; x end proc)); -1
dist := Distribution(Type = discrete, Support = 1 .. infinity, DiscreteValueMap = pontos, ProbabilityFunction = (proc (x) options operator, arrow; x end proc)); -1

> Z := RandomVariable(dist); -1

> map(convert, Sample(Z, 11), rational)
rtable(1 .. 11, [`/`(1, 4), `/`(1, 4), `/`(1, 2), `/`(1, 16), `/`(1, 4), `/`(1, 2), `/`(1, 2), `/`(1, 2), `/`(1, 4), `/`(1, 2), `/`(1, 8)], subtype = Vector[row]) (2.4.1.1.2)
 

> Mean(Z)
`/`(1, 3) (2.4.1.1.3)
 

> StandardDeviation(Z, numeric)
.1781741613 (2.4.1.1.4)
 

> plot([seq([pontos(n), pontos(n)], n = 1 .. 10)], style = point, symbolsize = 15)
Plot_2d
 

 

Definindo uma variável aleatória discreta a partir de uma amostra discreta com pesos 

Podemos definir uma variável aleatória discreta sobre uma amostra numérica da seguinte forma. Suponha que a amostra e os repectivos pesos sejam 

> Amostra := [1, 2, 3, 4]
[1, 2, 3, 4] (2.5.1)
 

> Pesos := [`/`(1, 2), `/`(1, 4), `/`(1, 8), `/`(1, 8)]
[`/`(1, 2), `/`(1, 4), `/`(1, 8), `/`(1, 8)] (2.5.2)
 

Para definir a variável aleatória discreta devemos usar o comando EmpliralDistribution dentro do comando RandomVariable da seguinte forma: 

Podemos agora usar qualquer comando do pacote Statistics que usem variáveis aleatórias. Por exemplo, amostragem, média e desvio padrão: 

> P := RandomVariable(EmpiricalDistribution(Amostra, probabilities = Pesos)); -1

> Sample(P, 10)

> Mean(P)

> StandardDeviation(P)
 

 

rtable(1 .. 10, [2., 2., 2., 2., 3., 2., 1., 1., 4., 2.], subtype = Vector[row])
(2.5.3)
 

 

Se tivemos uma amostra de números sem os pesos, podemos definir os pesos de cada elemento usando a frequência de ocorrência de cada elemento. Para definir uma variável aleatória baseada nestes pesos induzidos pela frequência, usamos também o comando EmpiricalDistribution da seguinte forma: 

> Amostra := [1, 1, 1, 1, 2, 2, 3, 4]
[1, 1, 1, 1, 2, 2, 3, 4] (2.5.4)
 

Podemos agora usar qualquer comando do pacote Statistics que usem variáveis aleatórias. Por exemplo, amostragem, média e desvio padrão: 

> P2 := RandomVariable(EmpiricalDistribution(Amostra)); -1

> Sample(P2, 10)

> Mean(P2)

> StandardDeviation(P2)
 

 

rtable(1 .. 10, [2., 1., 3., 3., 1., 4., 3., 2., 4., 2.], subtype = Vector[row])
(2.5.5)
 

Note que o segundo exemplo é exatamente igual ao primeiro. 

 

Exercício 

Suponha que jogamos dois dados. Seja X a variável aleatória dada pela soma dos valores dos dados. (a) Ache o valor médio e o desvio padrão de X. (b) Ache a distribuição de probabilidades de X e faça o histograma discreto da distribuição de probabilidades. 

Geração de números aleatórios sem o pacote Statistics 

O principal comando para gerar números inteiros aleatórios é o comando rand . Por exemplo, vamos gerar uma lista de números inteiros aleatórios: 

> randomize(); -1

> L := [seq((rand(100 .. 1000))(), i = 1 .. 30)]

(2.6.1)
 

A partir de L, podemos gerar uma lista de números racionais aleatórios entre `/`(1, 10) e 1, da seguinte forma: 

> L2 := map(proc (x) options operator, arrow; evalf(`+`(`*`(`/`(1, 1000), `*`(x))), 3) end proc, L)
[.801, .745, .835, .928, .895, .769, .851, .766, .881, .450, .810, .525, .938, .312, .367, .849, .851, .174, .859, .751, .578, .895, .989, .235, .618, .525, .510, .694, .331, .798]
[.801, .745, .835, .928, .895, .769, .851, .766, .881, .450, .810, .525, .938, .312, .367, .849, .851, .174, .859, .751, .578, .895, .989, .235, .618, .525, .510, .694, .331, .798]
(2.6.2)
 

A média e o desvio padrão podem ser calculados com os comando Mean e StandardDeviation do pacote Statistics: 

> Statistics[Mean](L2)
HFloat(0.6843333333333331) (2.6.3)
 

> Statistics[StandardDeviation](L2)
HFloat(0.22913199528736228) (2.6.4)
 

Um comando mais interessante para gerar estruturas com números aleatórios é Generate que está no pacote RandomTools . 

> with(RandomTools)
[AddFlavor, BlumBlumShub, Generate, GetFlavor, GetFlavors, GetState, HasFlavor, LinearCongruence, MersenneTwister, QuadraticCongruence, RemoveFlavor, SetState, returnValueInertProc]
[AddFlavor, BlumBlumShub, Generate, GetFlavor, GetFlavors, GetState, HasFlavor, LinearCongruence, MersenneTwister, QuadraticCongruence, RemoveFlavor, SetState, returnValueInertProc]
(2.6.5)
 

Para gerar uma lista (list) com 5 elementos do tipo ponto-flutuante (float) no intervalo (range) entre 10 a 20, o comando é 

> Generate(list(float(range = 10 .. 20), 5))
[15.45375030, 13.20173889, 13.20695677, 18.19148318, 12.44891140] (2.6.6)
 

Para gerar uma matrix `+`(`*`(3, `*`(x3))) com números aleatórios:  

> Matrix(3, 3, Generate(integer(range = 1 .. 10), makeproc = true))
rtable(1 .. 3, 1 .. 3, [[9, 8, 9], [10, 8, 10], [8, 5, 10]], subtype = Matrix) (2.6.7)
 

Média, Moda e Mediana 

Se tivermos uma sequência de números a[1] .. a[n], podendo ter repetições, a média ( Mean ) é definida como 

mu = `/`(`*`(sum(a[i], i = 1 .. n)), `*`(n)) 

A média também é usualmente chamada de expectância. 

 

A moda ( Mode ) é o número que mais aparece na sequência. Por exemplo, a moda de 1,2,2,5,6 é 2. Por outro lado, a sequência 1,3,3,5,5,6 tem duas modas: 3 e 5, neste caso é a sequência é chamada bimodal. Se houver mais do que duas modas, normalmente usamos o termo multimodal e se não há repetições, usamos o termo sequência amodal (O comando Mode retorna todos os elementos). Uma sequência que tem uma única moda também é chamada unimodal. 

 

A mediana ( Median ) é o valor numérico m tal que metade dos números da sequência são menores do que m e metade dos números são maiores que m. Por exemplo, a mediana de 1,2,2,5,6,7,8  é m = 5, pois três numeros (1,2,2) são menores do que 5 e três números (6,7,8) são maiores do que 5. Por outro lado, a mediana de 2,5,6,7 é m = 5.5, pois como m não pode estar na sequência ele é o valor médio entre 5 e 6.  

 

As definições de média, moda e mediana dadas acima usaram uma amostra de pontos, por isto, elas também são chamadas de média amostral, moda amostral e mediana amostral. Estes mesmos conceitos podem ser aplicados para variáveis aleatórias com distribuições de probabilidades discretas ou contínuas. Estas definições são estendidas devido ao contexto. Se X é uma variável aleatória em um espaço amostral discreto constituído pelo conjunto de pontos a[1] .. a[n] (sem repetições, n podendo ser infinito), a média E(X) é definida como 

E(X) = sum(`*`(P(X = a[i]), `*`(a[i])), i = 1 .. n) 

onde P(X = a[i]) é a probabilidade da variável X ter o valor a[i]. Se o espaço amostral é a reta real, a média é definida como 

E(X) = int(`*`(P(X = x), `*`(x)), x = `+`(`-`(infinity)) .. infinity) 

 

Quando o espaço amostral é a reta real, a moda é o número (ou números) na reta real cujo valor é um máximo global da distribuição de probabilidades. As vezes, máximos locais também são chamados de moda. 

 

A mediana é definida como o valor m tal que 

`>=`(P(`<=`(X, m)), `/`(1, 2)) e `>=`(P(`>=`(X, m)), `/`(1, 2)) 

O valor de m não é único necessariamente. Se `>`(P(X = x), 0) para todo x, o valor de m é único e geometricamente corresponde ao ponto da reta real que divide a área embaixo da curva da distribuição de probabilidades exatamente no meio.  

 

Por exemplo, a distribuição de probabilidades contínua `*`(`^`(chi, 2)) (chi-quadrado ou ChiSquare ) com nu = 5 graus de liberdade tem o seguinte gráfico 

>

> DensityPlot(ChiSquare(5), range = -3 .. 20, thickness = 2, color = blue)
DensityPlot(ChiSquare(5), range = -3 .. 20, thickness = 2, color = blue)
Plot_2d
 

Vamos calcular a média, moda e mediana para ChiSquare(5). 

> mu, moda, mediana := Mean(ChiSquare(5), numeric), Mode(ChiSquare(5), numeric), Median(ChiSquare(5), numeric)
mu, moda, mediana := Mean(ChiSquare(5), numeric), Mode(ChiSquare(5), numeric), Median(ChiSquare(5), numeric)
5., 3., HFloat(4.351461712458544) (2.7.1)
 

> x := 'x'; -1

Vamos conferir o resultado para a media usando a definição. 

> f := unapply(PDF(ChiSquare(5), x), x); -1

> int(`*`(f(x), `*`(x)), x = `+`(`-`(infinity)) .. infinity, numeric)
5.000000000 (2.7.2)
 

Vamos conferir o resultado para a moda. 

> evalf(maximize(f(x), x, location))
.1541803298, {[{x = 3.}, .1541803298]} (2.7.3)
 

Note que o máximo ocorre no ponto x = 3, que é a moda. Vamos conferir o resultado da mediana. 

> int(f(x), x = `+`(`-`(infinity)) .. mediana) = int(f(x), x = mediana .. infinity)
.5000002085 = .4999997915 (2.7.4)
 

Desvio Padrão, Variância, Assimetria e Curtose 

Se X é uma variável aleatória em um espaço amostral discreto constituído pelo conjunto de pontos a[1] .. a[n] (sem repetições, n podendo ser infinito) com média mu, a variância Var(X) ( Variance ) é definida como 

Var(X) = sum(`*`(P(X = a[i]), `*`(`^`(`+`(a[i], `-`(mu)), 2))), i = 1 .. n) 

onde P(X = a[i]) é a probabilidade da variável X ter o valor a[i]. Se o espaço amostral é a reta real, a média é definida como 

Var(X) = int(`*`(P(X = x), `*`(`^`(`+`(x, `-`(mu)), 2))), x = `+`(`-`(infinity)) .. infinity) 

Uma forma alternativa de definir variância é como a média E(Y) de uma nova variável aleatória Y definida como Y = `*`(`^`(`+`(X, `-`(mu)), 2)), isto é,  

Var(X) = E(`*`(`^`(`+`(X, `-`(mu)), 2))) 

Por exemplo, vamos checar que esta segunda definição usando uma variável aleatória X com a distribuição de probabilidades ChiSquare de 5 graus de liberdade: 

>

> X := RandomVariable(ChiSquare(5)); -1

> Variance(X) = Mean(`*`(`^`(`+`(X, `-`(Mean(X))), 2)))
10 = 10 (2.8.1)
 

O comando Variance quando aplicado em uma amostra ou lista de dados A[1], () .. (), A[n] calcula a variância usando estimativa imparcial ou não-viciada (unbiased estimation of variance), cuja fórmula é 

`/`(`*`(sum(`*`(`^`(`+`(A[i], `-`(conjugate(A))), 2)), i = 1 .. n)), `*`(`+`(n, `-`(1)))) 

 

O desvio padrão sigma(X) ( StandardDeviation ) é definido como 

sigma(X) = sqrt(Var(X)) 

O desvio padrão é uma medida do tamanho da dispersão do pontos de uma amostra em torno da média. Se a dispersão for pequena, o desvio padrão é pequeno e vice-versa. O comando StandardDeviation quando aplicado em uma amostra ou lista de dados usa a estimativa imparcial. Para encontrar a estimativa parcial, devemos multiplicar o resultado por sqrt(`/`(`*`(`+`(n, `-`(1))), `*`(n))). 

 

Exercícios 

1. Considere o gráfico de uma distribuição Normal de média mu e desvio padrão Mostre que a distância horizontal do ponto médio do gráfico até o ponto de inflexão (ponto onde a derivada segunda se anula) é sigma. 

 

2. É possível mostrar de maneira geral que  

E(`*`(`^`(`+`(X, `-`(mu)), 2))) = `+`(E(`*`(`^`(X, 2))), `-`(`*`(`^`(mu, 2)))) 

Verifique no Maple que este fato é verdeiro para a uma variável X com uma distribuição Normal de média mu e desvio padrão  

 

A assimetria ( Skewness ) é definida como 

alpha[3](X) = `/`(`*`(E(`*`(`^`(`+`(X, `-`(mu)), 3)))), `*`(`^`(sigma, 3))) 

A assimetria alpha[3]é uma medida da assimetria de uma distribuição de probabilidades: alpha[3] é positiva se for a distribuição de probabilidades for assimétrica à direita (veja o gráfico de `*`(`^`(chi(5), 2))), caso contrário alpha[3] é negativa. Por exemplo, a distribuição ChiSquare(nu) de nu graus de liberdade é assimétrica à equerda para qualquer valor de nu. Note que a assimetria alpha[3] é positiva. 

>

> nu := 'nu'; -1

> Skewness(ChiSquare(nu))
(2.8.2)
 

A distribuição Beta com o segundo parâmetro menor do que o primeiro é um exemplo de assimetria negativa 

> Skewness(BetaDistribution(5, 2))
`+`(`-`(`*`(`/`(4, 15), `*`(`^`(5, `/`(1, 2)))))) (2.8.3)
 

Note que o gráfico é assimétrico à esquerda: 

> DensityPlot(BetaDistribution(5, 2), range = -.2 .. 1.2)
Plot_2d
 

 

A curtose ( Kurtosis ) é definida como 

 

alpha[4](X) = `/`(`*`(E(`*`(`^`(`+`(X, `-`(mu)), 4)))), `*`(`^`(sigma, 4))) 

A curtose mede o grau de agudez do pico da distribuição de probabilidades tomando a distribuição Normal como referência. A curtose da distribuição Normal é 3. Se a curtose for maior do que 3, a distribuição é mais concentrada em torno da média do que a distribuição Normal e se for menor do que 3, a distribuiçãoé mais espalhada do que a distribuição Normal. Note que 

> mu, sigma := 'mu, sigma'; -1

> Kurtosis(Normal(mu, sigma))
3 (2.8.4)
 

O menor valor possível da curtose é 1 que é obtida com a distribuição de probabilidades discreta de Bernoulli quando a probabilidade de sucesso é `/`(1, 2). 

> Kurtosis(Bernoulli(`/`(1, 2)))
1 (2.8.5)
 

 

Exercício 

Moste usando Maple que a curtose da distribuição `*`(`^`(chi(nu), 2)) é sempre maior do que a curtose da distribuição normal. 

 

Percentil, Intervalo Interquartílico e Desvio Médio 

Considere uma distribuição de probabilidades contínua, por exemplo, a distribuição normal de média 0 e desvio padrão 1. Qual é o valor de x tal que a probabilidade de se obter um valor menor ou igual a x seja de 30%? A resposta seria o valor de x tal que a integral de `+`(`-`(infinity)) a x`*`(30, `/`(1, 100)). A interpretação geométrica é a área embaixo da curva da distribuição de probabilidades de `+`(`-`(infinity)) a x é `/`(3, 10). Este valor de x é chamado de trigésimo percentil. Para calcular um percentil, devemos usar o comando Percentile , no caso do trigésimo percentil, o comando é 

>

> x := 'x'; -1

> x[30] := Percentile(Normal(0, 1), 30, numeric); 1
(2.9.1)
 

Vamos verificar o resultado 

> t := 't'; -1

> int(PDF(Normal(0, 1), t), t = `+`(`-`(infinity)) .. x[30])

> x := 'x'; -1
(2.9.2)
 

Este valor de x também é chamado de terceiro decil, e também pode ser calculado com o comando Decile : 

> Decile(Normal(0, 1), 3, numeric); 1
HFloat(-0.5244005127080493) (2.9.3)
 

 

O intervalo interquartílico ( InterquartileRange ) é definido como a diferença entre o septagésimo-quinto (`^`(75, o)) percentil e o vigésimo-quinto (`^`(25, o)) percentil ou, equivalentemente, a diferença entre o terceiro quartil e o primeiro quartil ( Quartile ). Por exemplo, no caso da distribuição normal 

> InterquartileRange(Normal(0, 1), numeric)
HFloat(1.3489795003925509) (2.9.4)
 

Só para confirmar: 

> `+`(Percentile(Normal(0, 1), 75, numeric), `-`(Percentile(Normal(0, 1), 25, numeric)))
HFloat(1.3489795003925509) (2.9.5)
 

> `+`(Quartile(Normal(0, 1), 3, numeric), `-`(Quartile(Normal(0, 1), 1, numeric)))
HFloat(1.3489795003925509) (2.9.6)
 

O intervalo semi-interquartílico é definido como a metade do intervalo interquartílico. 

 

O desvio médio é definido como média E(Y) de uma nova variável aleatória Y definida como Y = abs(`+`(X, `-`(mu))), isto é,  

D(X) = E(abs(`+`(X, `-`(mu)))) 

> MeanDeviation(Normal(10, 1), numeric)
.7978845604 (2.9.7)
 

Podemos confirmar com os seguintes comandos: 

> X := RandomVariable(Normal(10, 1)); -1

> Y := abs(`+`(X, `-`(10))); -1

> Mean(Y, numeric)
(2.9.8)
 

 

O desvio mediano absoluto é definido como médiana de uma nova variável aleatória Y definida como Y = abs(`+`(X, `-`(mediana(X)))), isto é,  

MAD(X) = mediana(abs(`+`(X, `-`(mediana(X))))) 

Por exemplo: 

> MedianDeviation(Normal(0, 1), numeric)
HFloat(0.6744897501961059) (2.9.9)
 

Podemos confirmar com os seguintes comandos: 

> X := RandomVariable(Normal(10, 1)); -1

> Y := abs(`+`(X, `-`(Median(X)))); -1

> Median(Y, numeric)
(2.9.10)
 

Função Geradora de Momentos e Função Característica 

A partir de uma variável aleatória X podemos definir uma nova variável aleatória Y = f(X) usando uma função f. Podemos usar Y em os comandos que tem variáveis aleatórias como argumento, como o comando que calcula a média ou o desvio padrão.  Em particular, se a função for f(x) = exp(`*`(t, `*`(x))), teremos a nova variável aleatória Y defina como Y = exp(`*`(t, `*`(X))). A função geradora de momentos ( MomentGeneratingFunction ) é definida como a média de Y, isto é, 

M[X](t) = E(exp(`*`(t, `*`(X)))) 

O r-ésimo momento ( Moment ) de uma variável aleatória X e definido por  

(diff(mu(x), x))[r] = E(`^`(X, r)) 

e o r-ésimo momento central ( CentralMoment ) de uma variável aleatória X e definido como a média da variável aleatória `^`(`+`(X, `-`(mu)), r), onde mu é a média de X, isto é, por 

mu[r] = E(`^`(`+`(X, `-`(mu)), r)) 

Por exemplo 

>

> X := RandomVariable(Normal(10, 1)); -1

> F := unapply(MomentGeneratingFunction(X, t), t)
proc (t) options operator, arrow; exp(`+`(`*`(10, `*`(t)), `*`(`/`(1, 2), `*`(`^`(t, 2))))) end proc (2.10.1)
 

Note que se expandirmos a função geradora de momentos em série de Taylor 

> t := 't'; -1

> series(F(t), t)
series(`+`(1, `*`(10, `*`(t)), `*`(`/`(101, 2), `*`(`^`(t, 2))), `*`(`/`(515, 3), `*`(`^`(t, 3))), `*`(`/`(10603, 24), `*`(`^`(t, 4))), `*`(`/`(11015, 12), `*`(`^`(t, 5))))O(`^`(t, 6)),t,6) (2.10.2)
 

podemos obter o r-ésimo momento como o coeficiente do termo `^`(t, r) vezes o fatorial de r. De fato 

> mu := 'mu'; -1

> seq(mu[r] = `/`(`*`(Mean(`^`(X, r))), `*`(factorial(r))), r = 0 .. 5)
(2.10.3)
 

De outra maneira 

> for r to 5 do Moment(X, r) = ((`@@`(D, r))(F))(0) end do
 

 

 

 

10 = 10
101 = 101
1030 = 1030
10603 = 10603
110150 = 110150 (2.10.4)
 

> k := 11; -1

>

> X[i] := RandomVariable(EmpiricalDistribution([1, -1], probabilities = [`/`(1, 2), `/`(1, 2)])); -1

>

> expand(convert(MomentGeneratingFunction(add(X[i], i = 1 .. k), omega), trig))
(2.10.5)
 

Exercício 

Seja X[i] para i = 1 .. k variáveis aleatórias independentes que podem assumir os valores 1 ou -1 com probabilidade uniforme. Escolha k como sendo um número inteiro maior do que 10. (a) Mostre que a função geradora de momentos variável aleatória é `^`(cosh(t), k). (b) Mostre que a função geradora de momentos da variável aleatória é `^`(cosh(`/`(`*`(t), `*`(sqrt(k)))), k). 

Solução 

(a) Vamos mostrar no caso em que k = 11: 

>

> k := 11; -1

>

> X[i] := RandomVariable(EmpiricalDistribution([1, -1], probabilities = [`/`(1, 2), `/`(1, 2)])); -1

>

> expand(convert(MGF(add(X[i], i = 1 .. k), t), trig))
(2.10.1.1.1)
 

(b) 

> res := expand(convert(MGF(`/`(`*`(add(X[i], i = 1 .. k)), `*`(sqrt(k))), t), trig)); -1

> subs(T = `/`(`*`(t), `*`(k)), expand(subs(t = `*`(k, `*`(T)), res)))

> k := 'k'; -1
(2.10.1.1.2)
 

 

 

 

A função característica ( CharacteristicFunction ) é definida como a média de exp(`*`(I, `*`(omega, `*`(X)))), isto é, 

phi[X](omega) = E(exp(`*`(I, `*`(omega, `*`(X))))) 

Por exemplo 

> X := RandomVariable(Normal(10, 1)); -1

> F := unapply(CharacteristicFunction(X, omega), omega)
proc (omega) options operator, arrow; exp(`+`(`*`(`*`(10, `*`(I)), `*`(omega)), `-`(`*`(`/`(1, 2), `*`(`^`(omega, 2)))))) end proc (2.10.6)
 

Note que se expandirmos a função característica em série de Taylor 

> series(F(omega), omega)
series(`+`(1, `*`(`*`(10, `*`(I)), `*`(omega)), `-`(`*`(`/`(101, 2), `*`(`^`(omega, 2)))), `-`(`*`(`*`(`/`(515, 3), `*`(I)), `*`(`^`(omega, 3)))), `*`(`/`(10603, 24), `*`(`^`(omega, 4))), `*`(`*`(`/`(... (2.10.7)
 

podemos obter o r-ésimo momento como o coeficiente do termo `^`(t, r) vezes `/`(`*`(factorial(r)), `*`(`^`(I, r))). De fato 

> seq(mu[r] = `/`(`*`(Mean(`^`(`*`(I, `*`(X)), r))), `*`(factorial(r))), r = 0 .. 5)
mu[0] = 1, mu[1] = `*`(10, `*`(I)), mu[2] = -`/`(101, 2), mu[3] = `+`(`-`(`*`(`/`(515, 3), `*`(I)))), mu[4] = `/`(10603, 24), mu[5] = `*`(`/`(11015, 12), `*`(I)) (2.10.8)
 

De outra maneira 

> for r to 5 do Moment(X, r) = `/`(`*`(((`@@`(D, r))(F))(0)), `*`(`^`(I, r))) end do
 

 

 

 

10 = 10
101 = 101
1030 = 1030
10603 = 10603
110150 = 110150 (2.10.9)
 

Exercício 

Seja X[i] para i = 1 .. k variáveis aleatórias independentes que podem assumir os valores 1 ou -1 com probabilidade uniforme. Escolha k como sendo um número inteiro maior do que 10. (a) Mostre que a função característica da variável aleatória é `^`(cos(omega), k). (b) Mostre que a função característica da variável aleatória é `^`(cos(`/`(`*`(omega), `*`(sqrt(k)))), k). 

Estatística 

Amostragem 

Na Estatística, usualmente temos uma população de indivíduos (por exemplo: pessoas, produtos, ou objetos) e queremos tirar conclusões sobre as características desta população examinando um subgrupo que é chamado de amostra. Quando selecionamos uma amostra, vamos supor que cada indivíduo é selecionado aleatoriamente segundo uma distribuição de probabilidades associada à população. O tamanho da população pode ser grande, ou até mesmo infinito, mas o tamanho da amostra deve ser o menor possível por uma questão de economia de recursos. A amostragem pode ser com ou sem reposição, dependendo se um indivíduo selecionado foi colocado de volta na população ou não. No primeiro caso, amostra com reposição, um indivíduo pode ser selecionado mais do que uma vez. Os parâmetros da população são os parâmetros da distribuição de probabilidades associada. Por exemplo, uma população associada a distribuição normal (chamada população normal) tem como parâmetros a média mu e o desvio padrão sigma da distribuição normal. O objetivo da Estatística é usar dados amostrais para obter valores numéricos que sirvam para estimar e testar hipóteses sobre os parâmetros da população. Os valores numéricos são aproximações dos parâmetros da população e estimativas de erros são bem-vindos. 

 

Se fizermos várias amostragens, podemos definir novas variáveis aleatórias que terão suas próprias distribuições de probabilidades. Por exemplo, vamos supor que uma amostra X[1], () .. (), X[n] tem tamanho n. Vamos definir uma nova variável aleatória 

 

chamada média amostral. Podemos agora nos perguntar qual é a distribuição de probabilidades associada à Y uma vez conhecida a distribuição de probabilidades associada à X. Teremos um valor de Y cada vez que uma amostragem de tamanho n for realizada. Usuamente a variável Y é denotada por conjugate(X). 

 

Exemplo 1. A altura da população brasileira obedeçe aproximadamente uma distribuição normal de média mu = `+`(`*`(170, `*`(cm))) e desvio padrão sigma = `+`(`*`(7, `*`(cm))). Suponha que selecionamos amostras cada uma com 30 indivíduos. (a) Qual é a média e o desvio padrão da média amostral? (b) Suponha que selecionamos 50 amostras. Em quantas amostras encontramos a média entre `+`(`*`(160, `*`(cm))) e `+`(`*`(170, `*`(cm))). 

 

Solução. Antes de resolvermos este problema, note que podemos fazer uma simulação no Maple de um experimento deste tipo com a população brasileira.  

>

> Digits := 5

> UseHardwareFloats := false
 

(3.1.1)
 

> X := RandomVariable(Normal(170, 7)); -1

> for i from 0 to 50 do amostra[i] := Sample(X, 30) end do; -1

> L := [seq(Mean(amostra[i]), i = 1 .. 50)]
[170.17, 167.10, 171.27, 170.20, 172.24, 168.66, 168.10, 171.27, 172.22, 169.17, 168.69, 170.15, 169.91, 169.42, 170.05, 168.34, 169.91, 170.97, 168.49, 170.75, 168.14, 171.27, 171.48, 170.63, 169.54,...
[170.17, 167.10, 171.27, 170.20, 172.24, 168.66, 168.10, 171.27, 172.22, 169.17, 168.69, 170.15, 169.91, 169.42, 170.05, 168.34, 169.91, 170.97, 168.49, 170.75, 168.14, 171.27, 171.48, 170.63, 169.54,...
[170.17, 167.10, 171.27, 170.20, 172.24, 168.66, 168.10, 171.27, 172.22, 169.17, 168.69, 170.15, 169.91, 169.42, 170.05, 168.34, 169.91, 170.97, 168.49, 170.75, 168.14, 171.27, 171.48, 170.63, 169.54,...
[170.17, 167.10, 171.27, 170.20, 172.24, 168.66, 168.10, 171.27, 172.22, 169.17, 168.69, 170.15, 169.91, 169.42, 170.05, 168.34, 169.91, 170.97, 168.49, 170.75, 168.14, 171.27, 171.48, 170.63, 169.54,...
(3.1.2)
 

> Mean(L)
169.75 (3.1.3)
 

> StandardDeviation(L); 1
1.5080 (3.1.4)
 

> `*`(StandardDeviation(L), `*`(sqrt(`+`(`*`(24., `*`(`/`(25.)))))))
1.4775 (3.1.5)
 

Podemos ver que a média amostral deu muito próximo da média mu porém o desvio padrão das médias amostrais deu bem diferente de sigma. Os resultados corretos do item (a) podem ser obtidos da seguinte forma: 

> for i to 30 do X[i] := RandomVariable(Normal(170, 7)) end do; -1

> Y := `+`(`*`(`/`(1, 30), `*`(add(X[i], i = 1 .. 30)))); -1

> Mean(Y)
170 (3.1.6)
 

> evalf(StandardDeviation(Y))
1.2780 (3.1.7)
 

Pela teoria da Estatística, quando a população é infinita ou quando a amostragem é com reposição, o desvio padrão é dado por 

sigma[conjugate(X)] = `/`(`*`(sigma), `*`(sqrt(n))) 

No exemplo acima, temos `and`(sigma[conjugate(X)] = `+`(`/`(`*`(7), `*`(sqrt(30)))), `+`(`/`(`*`(7), `*`(sqrt(30)))) = 1.2781). Portanto, se fizermos um experimento com a população brasileira selecionando 50 amostras aleatórias de 30 indivíduos, após calcularmos sigma[conjugate(X)], podemos estimar o desvio padrão da população brasileira multiplicando o valor de sigma[conjugate(X)] dado pela Eq. evalf(StandardDeviation(Y)) por sqrt(n)=5.4772. 

 

O ítem (b) pede quantas amostras encontramos a média entre `+`(`*`(160, `*`(cm))) e `+`(`*`(170, `*`(cm))). Na simulação feita acima, o número de amostras é 

> nops(select(proc (x) options operator, arrow; `and`(`<=`(160, x), `<=`(x, 170)) end proc, L))
29 (3.1.8)
 

No entanto, o resultado correto é encontrado calculando a área embaixo da curva da distribuição de probabilidades de Y no intervalo [160, 170]. Isto é calculado da seguinte forma: 

> `+`(`*`(50, `*`(int(PDF(Y, x), x = 160 .. 170, numeric))))
25.000 (3.1.9)
 

> Digits := 10; -1

> UseHardwareFloats := true; -1

 

Um teorema importante da Estatística afirma que a variável aleatória padronizada Z = `/`(`*`(`+`(conjugate(X), `-`(mu[conjugate(X)]))), `*`(sigma[conjugate(X)])) é assintoticamente normal com média 0 e desvio padrão 1. Na prática, se o tamanho da amostra for `>=`(n, 30), o teorema pode ser aplicado com alguma margem de segurança. 

 

Exemplo 2. Uma máquina produz peças que pesam na média `+`(`*`(5, `*`(kg))) com um desvio padrão de `+`(`*`(.45, `*`(kg))). Um lote de 50 peças deve ser transportado. (a) Qual é a  probabilidade do peso total estar entre 245 e 255 kg? (b) Qual é a  probabilidade do peso total ser maior do que 255 kg?  

Solução. Vamos usar a distribuição normal com média mu = `+`(`*`(5, `*`(kg))) e desvio padrão sigma = `+`(`/`(`*`(.45, `*`(kg)), `*`(sqrt(n)))) onde n é o tamanho da amostra. A solução do ítem (a) é 

>

> mu := 5.0
5.0 (3.1.10)
 

> sigma := `+`(`/`(`*`(.45), `*`(sqrt(50.))))
0.6363961029e-1 (3.1.11)
 

> evalf(`*`(245, `/`(1, 50)) .. `*`(255, `/`(1, 50)))
4.900000000 .. 5.100000000 (3.1.12)
 

> int(PDF(Normal(mu, sigma), x), x = `*`(245, `/`(1, 50)) .. `*`(255, `/`(1, 50)))
.8838982561 (3.1.13)
 

e do ítem (b) é 

> int(PDF(Normal(mu, sigma), x), x = `*`(255, `/`(1, 50)) .. infinity)
0.5805087193e-1 (3.1.14)
 

 

Exercício 

O peso médio dos produtos de uma loja de departamentos é de `+`(`*`(15, `*`(kg))) e desvio padrão `+`(`*`(2.5, `*`(kg))).  Qual é a propabilidade de que 200 pacotes escolhidos aleatoriamente ultrapasse o peso limite do elevador, que é  

 

Suponha novamente que  X[1], () .. (), X[n] são variáveis aleatórias independentes associadas a mesma distribuição de probabilidades com média mu e desvio padrão sigma. A variância amostral  `*`(`^`(S, 2)) é definida como 

 

Pode-se mostrar que E(`*`(`^`(S, 2))) = `/`(`*`(`+`(n, `-`(1)), `*`(`^`(sigma, 2))), `*`(n)). Uma vez que `*`(`^`(S, 2)) é um estimador parcial da variância, definimos um novo estimador imparcial ou não-viciado 

`*`(`^`(\, 2)) = `/`(`*`(n, `*`(`^`(S, 2))), `*`(`+`(n, `-`(1))))  

que tem a propriedade E(`*`(`^`(\, 2))) = `*`(`^`(sigma, 2)). Portanto, `*`(`^`(\, 2)) é um estimador imparcial da variância. Note que os comandos Variance e StandardDeviation , quando aplicados em uma amostra, lista de dados ou matriz de dados, usam o estimamor imparcial para o cálculo da variância amostral e do desvio padrão amostral.  

 

Vamos verificar através de uma simulação para um valor pequeno de n que a média de `*`(`^`(S, 2)) é `/`(`*`(`+`(n, `-`(1)), `*`(`^`(sigma, 2))), `*`(n)).  

>

> n := 10; -1

> for i to n do X[i] := RandomVariable(ChiSquare(11)) end do; -1

> plot(PDF(ChiSquare(11), x), x = -1 .. 50)
Plot_2d
 

> Xbarra := `/`(`*`(add(X[i], i = 1 .. n)), `*`(n)); -1

> S2 := `/`(`*`(add(`*`(`^`(`+`(X[i], `-`(Xbarra)), 2)), i = 1 .. n)), `*`(n)); -1

> Mean(S2) = `/`(`*`(`+`(n, `-`(1)), `*`(`^`(StandardDeviation(X[2]), 2))), `*`(n))
`/`(99, 5) = `/`(99, 5) (3.1.15)
 

Para facilitar o cálculo da distribuição de probabilidades, um terceiro estimador é definido como 

 

A distribuição de probabilidades de `/`(`*`(n, `*`(`^`(S, 2))), `*`(`^`(sigma, 2))) é `*`(`^`(chi(`+`(n, `-`(1))), 2)), isto é, Chi-quadrado com `+`(n, `-`(1)) graus de liberdade, se as variáveis originais forem normais. Podemos verificar por simulações que a distribuição de probabilidades de `/`(`*`(n, `*`(`^`(S, 2))), `*`(`^`(sigma, 2))) é `*`(`^`(chi(`+`(n, `-`(1))), 2)) tomando valores específicos para n e definindo variáveis aleatórias normais. Por exemplo: 

>

> n := 5

> mu, sigma := 0, 1
 

(3.1.16)
 

> for i to n do X[i] := RandomVariable(Normal(mu, sigma)) end do; -1

> Xbarra := `/`(`*`(add(X[i], i = 1 .. n)), `*`(n)); -1

> nSsigma := `/`(`*`(add(`*`(`^`(`+`(X[i], `-`(Xbarra)), 2)), i = 1 .. n)), `*`(`^`(sigma, 2))); -1

> L := Sample(nSsigma, 500000); -1

> graf1 := Histogram(L, averageshifted = 5, color = blue); -1

> graf2 := DensityPlot(ChiSquare(`+`(n, `-`(1))), range = 0 .. `+`(`*`(3, `*`(n))), thickness = 3, color = red); -1

> display([graf1, graf2])
Plot_2d
 

Note que a curva vermelha, que é a distribuição de probabilidades de `*`(`^`(chi(4), 2)), é muito bem aproximada pelo histograma ( histogram ) de uma amostra com 50000 valores. Podemos testar facilmente para outros valores de n. 

 

Exemplo 3. Vimos no Exemplo 1 que a altura da população brasileira obedeçe aproximadamente uma distribuição normal de média mu = `+`(`*`(170, `*`(cm))) e desvio padrão sigma = `+`(`*`(7, `*`(cm))). Qual é o desvio padrão das variâncias usando o estimador imparcial para amostras de tamanho n = 10, isto é, sigma[`*`(`^`(\, 2))]?  

Solução.   

>

> n := 10; -1

> for i to n do X[i] := RandomVariable(Normal(170, 7)) end do; -1

> Xbarra := `/`(`*`(add(X[i], i = 1 .. n)), `*`(n)); -1

> Schapeu2 := `/`(`*`(add(`*`(`^`(`+`(X[i], `-`(Xbarra)), 2)), i = 1 .. n)), `*`(`+`(n, `-`(1)))); -1

> evalf(StandardDeviation(Schapeu2))
23.09882151 (3.1.17)
 

 

Exemplo 4. O desvio padrão dos pesos de uma população de uma cidade é sigma = `+`(`*`(15, `*`(kg))). Suponha que 200 pessoas são selecionadas aleatoriamente e o desvio padrão é calculado. Este processo é repetido diversas vezes. (a) Calcule a média da distribuição de desvios padrões. (b) Calcule o desvio padrão da distribuição de desvios padrões. 

 

Solução.  Vamos supor que a distribuição de probabilidades dos pesos da população desta cidade é normal. Portanto, A distribuição de probabilidades de `/`(`*`(n, `*`(`^`(S, 2))), `*`(`^`(sigma, 2))) é `*`(`^`(chi(`+`(n, `-`(1))), 2)). A variável aleatória S é dada por sqrt(`*`(`/`(`*`(`^`(nS, 2)), `*`(`^`(sigma, 2))), `/`(`*`(`^`(sigma, 2)), `*`(n)))). As respostas são encontradas calculando a média e o desvio padrão de S. 

(a) 

> n := 200; -1

> sigma := 15; -1

> nSsigma := RandomVariable(ChiSquare(`+`(n, `-`(1)))); -1

> S := sqrt(`/`(`*`(nSsigma, `*`(`^`(sigma, 2))), `*`(n))); -1

> evalf(Mean(S))
14.94366783 (3.1.18)
 

(b) 

> evalf(StandardDeviation(S))
.7495275792 (3.1.19)
 

Para amostras grandes (, o desvio padrão de S pode ser calculado através da fórmula sigma[S] = `/`(`*`(sigma), `*`(sqrt(`+`(`*`(2, `*`(n)))))). 

> evalf(`/`(`*`(sigma), `*`(sqrt(`+`(`*`(2, `*`(n)))))))
.7500000000 (3.1.20)
 

M. Spiegel, Probabilidade e Estatística, `^`(3, a) edição 

Cap 5, problemas 5.128 e 5.130. 

 

Estimação e intervalo de confiança 

Uma estatística é um estimador imparcial ou não viciado de um parâmetro se sua média é igual ao parâmetro correspondente na população. Se duas estimativas tem a mesma média, o estimador mais eficiente é aquele que tem o menor desvio padrão. Podemos ter estimativa por ponto, quando apenas um número é usado, ou estimativa por intervalo, quando mais de um número é usado. A confiança é a previsão de erro ou a precisão de uma estimativa.  

 

Suponha que uma estatística E  tem média mu[E]  e desvio padrão sigma[E]. Se a distribuição amostral de E  for aproximadamente normal, o valor de E  será encontrado no intervalo  

`+`(mu[E], `&+-`(`*`(z, `*`(sigma[E]))))  

com a probabilidade P(`and`(`<=`(`+`(`-`(`*`(z, `*`(sigma[E]))), mu[E]), E), `<=`(E, `+`(`*`(z, `*`(sigma[E])), mu[E])))), que é chamada de nível de confiança. O parâmetro z  é chamado valor crítico. Por exemplo, para uma distribuição normal (qualquer que seja a média e o desvio padrão), o nível de confiança é aproximadamente 95.45% quando o valor crítico é  z = 2.  

>

> mu, sigma := 'mu, sigma'
mu, sigma (3.2.1)
 

> z := 2
2 (3.2.2)
 

> evalf(int(PDF(Normal(mu, sigma), x), x = `+`(`-`(`*`(sigma, `*`(z))), mu) .. `+`(`*`(sigma, `*`(z)), mu)))
.9544997360 (3.2.3)
 

A área correspondente sob a curva da distribuição normal na variável padronizada é área azul do gráfico a seguir. 

> graf1 := plot(PDF(Normal(0, 1), x), x = `+`(`-`(z), `-`(3)) .. `+`(z, 3), color = blue, thickness = 3); -1

> graf2 := plot(PDF(Normal(0, 1), x), x = `+`(`-`(z)) .. z, filled = true, transparency = .8, color = blue); -1

> display([graf1, graf2])
Plot_2d
 

Usualmente o sub-índice de z denota a probabilidade acumulada até o valor z.  No caso acima, temos 

> evalf[4](CDF(Normal(0, 1), 2))
.9772 (3.2.4)
 

Assim, denotamos z[.9772] = 2. O valor de z  e do sub-índice para um nível de confiança de 95% é aproximadamente z[.975] = 1.96 e para um nível de confiança de 99% é aproximadamente z[.995] = 2.58. 

 

Intervalo de confiança para a média amostral para amostras grandes. 

Para amostras grandes (`>=`(n, 30)), a média amostal conjugate(X)  com um nível de confiança P tem o intervalo de confiança `+`(conjugate(X), `&+-`(`/`(`*`(z, `*`(sigma)), `*`(sqrt(n))))), onde sigma é o desvio padrão da população e n  é o tamanho da amostra. Se o desvio padrão da população não é conhecido, usamos o intervalo `+`(conjugate(X), `&+-`(`/`(`*`(z, `*`(\)), `*`(sqrt(n))))), onde \  é o estimador imparcial do desvio padrão. O valor crítico z é calculado usando a distribuição normal, pois para n  grande a variável aleatória padronizada Z = `/`(`*`(`+`(conjugate(X), `-`(mu[conjugate(X)]))), `*`(sigma[conjugate(X)])) é assintoticamente normal com média 0 e desvio padrão 1. 

 

Continuação do Exemplo 1. A altura da população brasileira obedeçe aproximadamente uma distribuição normal de média mu = `+`(`*`(170, `*`(cm))) e desvio padrão sigma = `+`(`*`(7, `*`(cm))). Suponha que selecionamos 50 amostras cada uma com 30 indivíduos. (a) Qual é a média e o desvio padrão das 50 amostras? (b) Em quantas amostras encontramos a média entre `+`(`*`(160, `*`(cm))) e `+`(`*`(170, `*`(cm))). (c) Qual é o intervalo para se ter 95% de confiança para estimar a altura média. 

 

Solução. Os ítens (a) e (b) foram resolvidos na seção anterior. Vamos agora resolver o ítem (c). Quando o desvio padrão é conhecido, o intervalo de confiança é dado por `+`(conjugate(X), `&+-`(`/`(`*`(z, `*`(sigma)), `*`(sqrt(n))))). 

>

> z := fsolve(CDF(Normal(0, 1), x) = .975)
1.959963985 (3.2.5)
 

Usamos o valor 0.975 pois devemos acrescentar a área de 2.5% aos 95% para poder usar a função de distribuição acumulada. 

> sigma, n := 7., 30
7., 30 (3.2.6)
 

> mu_a, mu_b := evalf(`+`(170., `-`(`/`(`*`(z, `*`(sigma)), `*`(sqrt(n)))))), evalf(`+`(170., `/`(`*`(z, `*`(sigma)), `*`(sqrt(n)))))
167.4951282, 172.5048718 (3.2.7)
 

> `and`(`<=`(mu_a, mu), `<=`(mu, mu_b))
`and`(`<=`(167.4951282, mu), `<=`(mu, 172.5048718)) (3.2.8)
 

Uma maneira de fazer uma simulação para testar este intervalo é 

> N := 5000; -1

> X := RandomVariable(Normal(170, 7)); -1

> for i to N do amostra[i] := Sample(X, n) end do; -1

> L := [seq(Mean(amostra[i]), i = 1 .. N)]; -1

> `*`(evalf(`+`(`/`(`*`(100), `*`(N)))), `*`(nops(select(proc (mu) options operator, arrow; `and`(`<=`(mu_a, mu), `<=`(mu, mu_b)) end proc, L))))
(3.2.9)
 

Note que o resultado acima tem que dar próximo de 95%. Quanto maior for o número de amostras, mais perto de 95%. 

 

 

Exercício 

Calcule os valores críticos para um nível de confiança de 95% para a distribuição t de Student. 

 

 

Intervalo de confiança para a média amostral para amostras pequenas. 

Para amostras pequenas (`<`(n, 30)), a média amostal conjugate(X)  com um nível de confiança P tem o intervalo de confiança `+`(conjugate(X), `&+-`(`/`(`*`(t, `*`(sigma)), `*`(sqrt(n))))), onde sigma é o desvio padrão de população normal e n  é o tamanho da amostra. Se o desvio padrão da população normal não é conhecido, usamos o intervalo `+`(conjugate(X), `&+-`(`/`(`*`(t, `*`(\)), `*`(sqrt(n))))), onde \  é o estimador imparcial do desvio padrão. O valor crítico t é calculado usando a distribuição t de Student (StudentT ), pois para n  pequeno a variável aleatória padronizada Z = `/`(`*`(`+`(conjugate(X), `-`(mu[conjugate(X)])), `*`(sqrt(n))), `*`(\)) obedece a distribuição t de Student. Portanto, para um nível de confiança de 95%, a variável Z  deve estar no invervalo `and`(`<`(`+`(`-`(t[.975])), Z), `<`(Z, t[.975])), que produz o intervalo `+`(conjugate(X), `&+-`(`/`(`*`(t, `*`(\)), `*`(sqrt(n))))). 

 

Vamos fazer uma simulação para checar que a variável aleatória padronizada Z = `/`(`*`(`+`(conjugate(X), `-`(mu[conjugate(X)])), `*`(sqrt(n))), `*`(\)) obedece a distribuição t de Student e comparar com a distribuição normal. Vamos considerar amostras de tamanho n = 10. 

>

> n, mu := 10, 0; -1

> for i to n do X[i] := RandomVariable(Normal(mu, 1)) end do; -1

> Xbarra := `/`(`*`(add(X[i], i = 1 .. n)), `*`(n)); -1

> S2 := `/`(`*`(add(`*`(`^`(`+`(X[i], `-`(Xbarra)), 2)), i = 1 .. n)), `*`(`+`(n, `-`(1)))); -1

> T := `/`(`*`(`+`(Xbarra, `-`(mu)), `*`(sqrt(n))), `*`(sqrt(S2))); -1

> L := Sample(T, 100000); -1

> graf1 := Histogram(L, averageshifted = 3, color = blue); -1

> distT := Distribution(StudentT(`+`(n, `-`(1)))); -1

> graf2 := DensityPlot(distT, range = -5 .. 5, color = red, legend = [

> graf3 := DensityPlot(Normal(0, StandardDeviation(L)), range = -5 .. 5, color = green, legend = [
graf3 := DensityPlot(Normal(0, StandardDeviation(L)), range = -5 .. 5, color = green, legend = [

> display([graf1, graf2, graf3])
Plot_2d
 

Para um nível de confiança de 95%, o valor de t[.975]  é  

> t := 't'; -1

> t[.975] := fsolve(CDF(StudentT(9), x) = .975)
(3.2.10)
 

que corresponde ao valor no eixo horizontal que delimita a área hachuriada no gráfico abaixo. 

> graf1 := plot(PDF(StudentT(9), t), t = -5 .. 5, color = blue, thickness = 3, tickmarks = [[t[.975] =
graf1 := plot(PDF(StudentT(9), t), t = -5 .. 5, color = blue, thickness = 3, tickmarks = [[t[.975] =

> graf2 := plot(PDF(StudentT(9), t), t = -5 .. t[.975], filled = true, transparency = .8, color = blue); -1

> display([graf1, graf2])
Plot_2d
 

 

Exercício 

Uma pessoa mediu o peso de uma peça em 12 balanças diferentes obtendo uma média de `+`(`*`(5.67, `*`(kg))) e desvio padrão de `+`(`*`(0.9e-1, `*`(kg))). Encontre os limites de 95% de confiança para o peso real da peça. 

Solução 

>

> t[.975] := fsolve(CDF(StudentT(11), x) = .975)
2.200985160 (3.2.2.1.1)
 

O peso real da peça está no intervalo 

> evalf([`+`(5.67, `-`(`/`(`*`(0.9e-1, `*`(t[.975])), `*`(sqrt(11))))), `+`(5.67, `/`(`*`(0.9e-1, `*`(t[.975])), `*`(sqrt(11))))])
[5.610274020, 5.729725980] (3.2.2.1.2)
 

com um nível de confiança de 95%. Usando o pacote Tolerances ,  

podemos expressar o intervalo da seguinte forma: 

> with(Tolerances); -1

> intervalo := `&+-`(5.67, `+`(`/`(`*`(0.9e-1, `*`(t[.975])), `*`(sqrt(11)))))
INTERVAL(5.610274020 .. 5.729725980) (3.2.2.1.3)
 

 

 

Intervalo de confiança para a variância de uma distribuição normal. 

Uma vez que `/`(`*`(n, `*`(`^`(S, 2))), `*`(`^`(sigma, 2))) = `/`(`*`(`+`(n, `-`(1)), `*`(`^`(\, 2))), `*`(`^`(sigma, 2))) tem um distribuição `*`(`^`(chi(`+`(n, `-`(1))), 2)), os limites de confiança para sigma é obtido a partir do intervalo 

`and`(`<=`(`*`(`^`(chi, 2), `*`([0.25e-1])), `/`(`*`(`+`(n, `-`(1)), `*`(`^`(\, 2))), `*`(`^`(sigma, 2)))), `<=`(`/`(`*`(`+`(n, `-`(1)), `*`(`^`(\, 2))), `*`(`^`(sigma, 2))), `*`(`^`(chi, 2), `*`([.97... 

para um nível de confiança de 95%. Usando o intervalo acima, obtemos 

`and`(`<=`(`/`(`*`(\, `*`(sqrt(`+`(n, `-`(1))))), `*`(`^`(chi, 2), `*`([.975]))), sigma), `<=`(sigma, `/`(`*`(\, `*`(sqrt(`+`(n, `-`(1))))), `*`(`^`(chi, 2), `*`([0.25e-1]))))) 

Os valores de `*`(`^`(chi, 2), `*`([0.25e-1]))  e `*`(`^`(chi, 2), `*`([.975]))  para 10 graus de liberdade são 

>

> fsolve(CDF(ChiSquare(10), x) = 0.25e-1)
3.246972780 (3.2.11)
 

> fsolve(CDF(ChiSquare(10), x) = .975)
20.48317735 (3.2.12)
 

 

Exercício 

Um lote de 150 lâmpadas foi analisado e o desvio padrão dos tempos de vida útil foi 80 horas. Encontre o limite de confiaça de 95% para o desvio padrão. 

Solução 

> z[.975] := fsolve(CDF(Normal(0, 1), x) = .975)
1.959963985 (3.2.3.1.1)
 

O desvio padrão está no intervalo 

> evalf([`+`(80, `-`(`/`(`*`(80, `*`(z[.975])), `*`(sqrt(`*`(2, 150)))))), `+`(80, `/`(`*`(80, `*`(z[.975])), `*`(sqrt(`*`(2, 150)))))])
[70.94731412, 89.05268588] (3.2.3.1.2)
 

com um nível de confiança de 95%. Usando o pacote Tolerances ,  

podemos expressar o intervalo da seguinte forma: 

> with(Tolerances); -1

> intervalo := `&+-`(80, `+`(`/`(`*`(80, `*`(z[.975])), `*`(sqrt(`*`(2, 150))))))
INTERVAL(70.94731412 .. 89.05268588) (3.2.3.1.3)
 

Note que o comando &+- pode ser usado fora do pacote Tolerances . Outro detalhe é que &+- não precisa estar entre crases, como é usual quando usamos caracteres não alfanuméricos em nomes. O caracter ampersand (&) é usado pelo usuário na definição de novos operadores no notação infix. Estes operadores também podem ser usados na notação prefix, porém neste caso as crases são necessárias, por exemplo 

> `&+-`(10, 0.1e-1)
INTERVAL(9.99 .. 10.01) (3.2.3.1.4)
 

 

Teste de Hipóteses e Significância 

Para tomar decisões baseadas em dados estatísticos, formulamos uma hipótese cujo grau de significância deve testado. Quando formulamos uma hipótese sobre um valor de uma estatística, ela é chamada de hipótese nula quando a hipótese afirma que o valor é o valor médio da estatística. As hipóteses que afirmam que é a média é diferente (maior ou menor) que o valor escolhido na hipótese nula são chamadas de hipóteses alternativas. Se rejeitamos uma hipótese nula verdadeira, dizemos que cometemos um erro do tipo I (falso positivo). Se aceitarmos uma hipótese nula falsa, dizemos que cometemos um erro do tipo II (falso negativo). Para testar uma hipótese, a probabilidade máxima com a qual estamos arriscando cometer um erro do tipo I é chamada de nível de significância. Na prática, usamos um nível de significância de 5%. Para questões sensíveis, usamos um nível de significância de 1%.  

 

Exemplo 1. Temos uma moeda queremos testar a seguinte hipótese nula: a moeda é não-viciada. O teste consiste em jogar a moeda 100 vezes e anotar quantas vezes deu cara ou coroa. Como regra de decisão vamos seguir a seguinte estratégia: (1) aceitar a hipótese se o número de caras está entre 40 e 60 inclusive, (2) rejeitar a hipótese, caso contrário. Encontre o nível de significância desta regra de decisão? 

Solução: Vamos calcular a probabilidade de rejeitar a hipótese quando ela está correta (erro do tipo I). A probabilidade de se conseguir entre 40 e 60 caras incluindo os extremos em 100 lançamento da moeda é 

> p1 := add(ProbabilityFunction(Binomial(100, .5), i), i = 40 .. 60)
.9647997995 (3.3.1)
 

A probabilidade de se cometer um erro do tipo I é 

> prob_erro_tipo_I := `+`(1, `-`(p1))
0.352002005e-1 (3.3.2)
 

O nível de significância é portanto aproximadamente 3.52%. 

 

Métodos Bayesianos 

Teorema de Bayes, probabilidades a priori e a posteriori 

 

Suponha que A é um evento e seja conjugate(A) o seu complemento. Considere um novo evento B. Então, o teorema de Bayes afirma que  

 

 

 

onde é a probabilidade condicional de ocorrer o evento A dado que o evento B já ocorreu. 

 

Exemplo 1. Suponha que temos duas moedas honestas (H) e uma moeda viciada (V), cuja probabilidade de dar cara é `/`(1, 5). Uma pessoa pega uma delas ao acaso e a joga três vezes. Suponha que o resultado divulgado foi duas caras e uma coroa (ca2co1) sem informação sobre a ordem. Qual é a probabilidade da moeda selecionada ser honesta? 

Solução. Pelo teorema de Bayes temos que 

 

onde  é a probabilidade da moeda escolhida ser honesta condicionado ao resultado ter sido duas caras e uma coroa (o evento ca2co1 já ocorreu). Portanto  é a resposta do problema. Temos que calcular cada probabilidade que aparece no lado direito da equação. Vamos começar pelas mais simples.  

 

> `P(H)` := `/`(2, 3); 1

> `P(V)` := `/`(1, 3); 1

> `P(ca2co1|H)` := `+`(`/`(`*`(3), `*`(`^`(2, 3)))); 1
`P(ca2co1|H)` := `+`(`/`(`*`(3), `*`(`^`(2, 3)))); 1

> `P(ca2co1|V)` := `+`(`*`(`/`(4, 5), `*`(`*`(`/`(3, 5), `/`(1, 5))))); 1
`P(ca2co1|V)` := `+`(`*`(`/`(4, 5), `*`(`*`(`/`(3, 5), `/`(1, 5))))); 1
 

 

 

(3.4.1.1)
 

 

A resposta do problema é 

> `P(H|ca2co1)` = evalf(`/`(`*`(`P(H)`, `*`(`P(ca2co1|H)`)), `*`(`+`(`*`(`P(H)`, `*`(`P(ca2co1|H)`)), `*`(`P(V)`, `*`(`P(ca2co1|V)`))))))
`P(H|ca2co1)` = .8865248227 (3.4.1.2)
 

Após saber o resultado (2 caras e 1 coroa), ficamos mais confiantes de que a moeda é honesta. Sem saber do resultado só poderíamos dizer que a probabilidade da moeda selecionada ser honesta era P(H) = .667. A probabilidade P(H) é chamada probabilidade a priori e é obtida antes de quaisquer testes ou observações. A probabilidade é chamada probabilidade a posteriori, pois é calculada após testes e observações terem sido feitos. 

 

Exemplo 2. Suponha que temos várias moedas honestas (H) e várias moedas viciadas (V), cuja probabilidade de dar cara é `/`(1, 5). Uma pessoa pega uma delas ao acaso e a joga três vezes. Suponha que o resultado divulgado foi duas caras e uma coroa (ca2co1) sem informação sobre a ordem. Qual é a probabilidade da moeda selecionada ser honesta? 

Solução. Não temos como resolver este problema pelo teorema de Bayes, pois não conhecemos nem P(H) nem P(V). O método Bayesiano supõe que o evento theta (que pode ser H ou V) antes de qualquer teste obedece à uma distribuição de probabilidades Pi(theta). Quando não há nenhuma informação disponível, a distribuição constante é usada, neste caso Pi(H) = `/`(1, 2) e Pi(V) = `/`(1, 2). Após aplicar o teorema da Bayes 

 

 

com `and`(P(H) = P(V), P(V) = `/`(1, 2)), temos 

> `P(H)` := `/`(1, 2); 1

> `P(V)` := `/`(1, 2); 1

> `P(ca2co1|H)` := `+`(`/`(`*`(3), `*`(`^`(2, 3)))); 1
`P(ca2co1|H)` := `+`(`/`(`*`(3), `*`(`^`(2, 3)))); 1

> `P(ca2co1|V)` := `+`(`*`(`/`(4, 5), `*`(`*`(`/`(3, 5), `/`(1, 5))))); 1
`P(ca2co1|V)` := `+`(`*`(`/`(4, 5), `*`(`*`(`/`(3, 5), `/`(1, 5))))); 1

> `P(H|ca2co1)` = eval(`/`(`*`(`P(H)`, `*`(`P(ca2co1|H)`)), `*`(`+`(`*`(`P(H)`, `*`(`P(ca2co1|H)`)), `*`(`P(V)`, `*`(`P(ca2co1|V)`))))))
 

 

 

 

(3.4.1.3)
 

Após saber o resultado (2 caras e 1 coroa), podemos atualizar nosso conhecimento sobre a distribuição de probabilidades Pi(theta). Agora usamos os valores  Pi(H) = `/`(125, 157) e Pi(V) = `/`(32, 157).  

 

 

O teorema (ou regra) de Bayes na sua forma mais geral é formulado da seguinte maneira. Suponha que A[1], () .. (), A[n] sejam eventos mutuamente exclusivos cuja união cobre o espaço amostral, isto é, sum(P(A[i]), i = 1 .. n) = 1. Considere um novo evento B no mesmo espaço amostral. Então, o teorema de Bayes afirma que  

 

 

 

 

Se o espaço amostral for contínuo, por exemplo, a reta real, a versão contínua do teorema de Bayes é 

 

 

onde theta[0] é um evento fixo do espaço amostral (um ponto na reta), é a probabilidade condicional de ocorrer theta[0] dado que x ocorreu e é a probabilidade condicional de ocorrer x dado que theta ocorreu A probabilidade a priori de ocorrer o evento theta[0] é Pi(theta[0]) e a probabilidade a posteriori é condicionado ao resultado ter sido x. 

 

Exemplo 3 (espaço amostral discreto). Suponha que a variável aleatória X tem a distribuição normal de média 0 e desvio padrão sigma desconhecido, porém os possíveis valores são sigma = 1 com probabilidade a priori Pi(1) = `/`(1, 2), sigma = 2 com probabilidade a priori Pi(2) = `/`(1, 4), sigma = 3 com probabilidade a priori Pi(3) = `/`(1, 8), e sigma = 4 com probabilidade a priori Pi(4) = `/`(1, 8). Uma amostra aleatória de X de tamanho 3 retornou [2.5, 3.1, 0.7]. Determine a distribuição a posteriori. 

Solução. As probabilidades a priori são 

>

> pi[1], pi[2], pi[3], pi[4] := `/`(1, 2), `/`(1, 4), `/`(1, 8), `/`(1, 8)
`/`(1, 2), `/`(1, 4), `/`(1, 8), `/`(1, 8) (3.4.1.4)
 

Suponha que  X tem a distribuição normal de média 0 e desvio padrão sigma. Portanto a probabilidade condicional da amostra ter sido [2.5, 3.1, 0.7] condicionado ao desvio padrão ser sigma é  

> sigma := 'sigma'; -1

> fxsigma := unapply(`*`(PDF(Normal(0, sigma), 2.5), `*`(PDF(Normal(0, sigma), 3.1), `*`(PDF(Normal(0, sigma), .7)))), sigma)
fxsigma := unapply(`*`(PDF(Normal(0, sigma), 2.5), `*`(PDF(Normal(0, sigma), 3.1), `*`(PDF(Normal(0, sigma), .7)))), sigma)
(3.4.1.5)
 

Usando o teorema de Bayes discreto, temos que a probabilidade a posteriori é calculada através da seguinte função 

> p_a_posteriori := unapply(`/`(`*`(pi[sigma0], `*`(fxsigma(sigma0))), `*`(add(`*`(pi[sigma], `*`(fxsigma(sigma))), sigma = 1 .. 4))), sigma0)
proc (sigma0) options operator, arrow; `+`(`/`(`*`(138.3626935, `*`(pi[sigma0], `*`(exp(`+`(`-`(`/`(`*`(3.125000000), `*`(`^`(sigma0, 2)))))), `*`(exp(`+`(`-`(`/`(`*`(4.805000000), `*`(`^`(sigma0, 2))... (3.4.1.6)
 

A probabilidades a posteriori correspondentes ao valores sigma = 1, sigma = 2, sigma = 3, sigma = 4 são 

> L := [seq(p_a_posteriori(sigma), sigma = 1 .. 4)]
[0.1948191402e-1, .5601181973, .2582742460, .1621256425] (3.4.1.7)
 

Vamos confirmar que L é uma distribuição de probabilidades discretas 

> add(i, i = L)
.9999999998 (3.4.1.8)
 

 

 

Exemplo 4 (espaço amostral contínuo). Suponha que a variável aleatória X tem a distribuição Normal de média theta desconhecida e desvio padrão sigma = 1. Suponha o valor de theta é também aleatório cuja probabilidade a priori é dada pela distribuição Normal padronizada (média 0 e desvio padrão 1).  Foi feita uma amostragem de X de tamanho 1 cujo resultado foi `/`(1, 2). Qual é a média da distribuição a posteriori? 

Solução.  A variável aleatória X tem a distribuição Normal de média theta e desvio padrão sigma = 1. Seja f  a distribuição de probabilidades de X que depende do valor de theta, assim 

>

> f := unapply(PDF(Normal(theta, 1), `/`(1, 2)), theta); 1
proc (theta) options operator, arrow; `+`(`/`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(`+`(`/`(1, 2), `-`(theta)), 2))))))))), `*`(`^`(Pi, `/`(1, 2))))) end proc (3.4.1.9)
 

A distribuição a priori (a distribuição de probabilidades de theta) é a distribuição Normal de média 0 e desvio padrão 1.  

Pelo teorema de Bayes (versão contínua), a distribuição a posteriori é 

> pi_a_priori := proc (theta) options operator, arrow; PDF(Normal(0, 1), theta) end proc; -1

> pi_a_posteriori := unapply(simplify(`/`(`*`(pi_a_priori(theta0), `*`(f(theta0))), `*`(int(`*`(pi_a_priori(theta), `*`(f(theta))), theta = `+`(`-`(infinity)) .. infinity)))), theta0)
proc (theta0) options operator, arrow; `/`(`*`(exp(`+`(`-`(`*`(`/`(1, 16), `*`(`^`(`+`(`*`(4, `*`(theta0)), `-`(1)), 2))))))), `*`(`^`(Pi, `/`(1, 2)))) end proc (3.4.1.10)
 

A média da distribuição a posteriori é 

> Mean(RandomVariable(Distribution(PDF = pi_a_posteriori)))
`/`(1, 4) (3.4.1.11)
 

 

 

Exemplo 5. Suponha que foi realizada uma pesquisa de opinião pública para determinar se o resultado do plebliscito seria SIM ou NÃO. De 1000 pessoas entrevistadas, 480 optarão pelo SIM e 520 pelo NÃO. Quais são as chances do SIM vencer? 

Solução. Uma vez que só há duas possibilidades, seja X uma variável aleatória com a distribuição binomial, cuja expressão analítica é 

>

> n, x := 'n, x'; -1

> expr := `assuming`([ProbabilityFunction(Binomial(n, theta), x)], [`>`(x, 0)])
(3.4.1.12)
 

onde n é o tamanho da amostra, x é o número de sucessos e theta é a probabilidade de sucesso (SIM). Queremos deteminar o valor de Para resolver este problema usando o método Bayesiano, vamos supor inicialmente que theta tem uma distribuição a priori uniforme no intervalo real [0,1],  

> pi := proc (theta) options operator, arrow; piecewise(`and`(`<=`(0, theta), `<`(theta, 1)), 1, 0) end proc
proc (theta) options operator, arrow; piecewise(`and`(`<=`(0, theta), `<`(theta, 1)), 1, 0) end proc (3.4.1.13)
 

Isto reflete nosso desconhecimento completo sobre a preferência do público antes da pesquisa de opinião. Queremos encontrar a distribuição a posteriori. Antes de usar o teorema de Bayes, temos que encontrar . Como X tem a distribuição binomial, n = 1000, temos que 

> f := unapply(expr, x, theta)
proc (x, theta) options operator, arrow; `*`(binomial(n, x), `*`(`^`(theta, x), `*`(`^`(`+`(1, `-`(theta)), `+`(n, `-`(x)))))) end proc (3.4.1.14)
 

Usando o teorema de Bayes na versão contínua 

> prob_posteriori := unapply(`/`(`*`(pi(theta0), `*`(f(x, theta0))), `*`(int(`*`(pi(theta), `*`(f(x, theta))), theta = 0 .. infinity))), theta0, x)
proc (theta0, x) options operator, arrow; `/`(`*`(piecewise(`and`(`<=`(0, theta0), `<`(theta0, 1)), 1, 0), `*`(`^`(theta0, x), `*`(`^`(`+`(1, `-`(theta0)), `+`(n, `-`(x))), `*`(GAMMA(`+`(n, 2)))))), `... (3.4.1.15)
 

O gráfico em função de theta quando a amostra tem tamanho 1000 e 480 sucessos é 

> n, x := 1000, 480
1000, 480 (3.4.1.16)
 

> plot(evalf(prob_posteriori(theta, x)), theta = 0 .. 1)
Plot_2d
 

O SIM vence se `and`(`<`(`/`(1, 2), theta), `<=`(theta, 1)). A probabilidade disto acontecer é 

> `+`(`*`(100, `*`(evalf(int(prob_posteriori(theta, x), theta = `/`(1, 2) .. 1)), `*`(`% `))))
`+`(`*`(10.30536554, `*`(`% `))) (3.4.1.17)
 

 

Vamos tentar reconhecer a expressão da probabilidade a posteriori. A expressão analítica é 

> evalf(prob_posteriori(theta, x))
`+`(`/`(`*`(0.5318282738e-1080, `*`(piecewise(`and`(`<=`(0., theta), `<`(theta, 1.)), 1., 0.), `*`(`^`(theta, 480), `*`(`^`(`+`(1., `-`(`*`(1., `*`(theta)))), `+`(n, `-`(480.))), `*`(GAMMA(`+`(n, 2.))... (3.4.1.18)
 

Compare este resultado com a expressão analítica da distribuição contínua Beta com parametros a e b. 

> convert(PDF(Beta(a, b), theta), GAMMA)
piecewise(`<`(theta, 0), 0, `<`(theta, 1), `/`(`*`(`^`(theta, `+`(`-`(1), a)), `*`(`^`(`+`(1, `-`(theta)), `+`(`-`(1), b)), `*`(GAMMA(`+`(a, b))))), `*`(GAMMA(a), `*`(GAMMA(b)))), 0) (3.4.1.19)
 

Facilmente determinamos os valores dos parametros a e b. 

> simplify(`+`(prob_posteriori(theta, x), `-`(convert(subs(a = `+`(x, 1), b = `+`(n, `-`(x), 1), PDF(Beta(a, b), theta)), GAMMA))))
0 (3.4.1.20)
 

Portanto, a distribuição a posteriori é a distribuição Beta com parâmetros `+`(x, 1) e `+`(n, `-`(x), 1), isto é, número de SIM mais 1 e número de NÃO mais 1. 

 

 

A distribuição de Poisson (pronúncia: pôa-som) é uma distribuição de probabilidades discreta para n inteiro não-negativo que depende de um parâmetro real positivo lambda. A expressão analítica desta distribuição é 

> P := `assuming`([evalf(ProbabilityFunction(Poisson(lambda), n))], [`>`(n, 0)])
`+`(`*`(0.2485168143e-2567, `*`(`^`(lambda, 1000), `*`(exp(`+`(`-`(`*`(1., `*`(lambda))))))))) (3.4.1.21)
 

O gráfico da distribuição são os pontos vermelhos. 

> n := 'n'; -1

> PoissonGraph := proc (l) options operator, arrow; display([plot([subs(lambda = l, P)], n = 0 .. 20, color = blue), plot([seq([i, evalf(subs(n = i, lambda = l, P))], i = 1 .. 20)], style = point, symbo...
PoissonGraph := proc (l) options operator, arrow; display([plot([subs(lambda = l, P)], n = 0 .. 20, color = blue), plot([seq([i, evalf(subs(n = i, lambda = l, P))], i = 1 .. 20)], style = point, symbo...

> display([PoissonGraph(1), PoissonGraph(4), PoissonGraph(10)]); 1
Plot_2d
 

A média de uma distribuição de Poisson de parâmetro lambda é exatamente lambda. Podemos confirmar com o seguinte comando. 

> Mean(Poisson(lambda))
lambda (3.4.1.22)
 

A distribuição de Poisson pode ser usada em muitas situações. Um exemplo típico envolve número de acidentes por ano em estradas. O número de acidente em 2014 foi de aproximadamente seja 169 mil. Podemos modelar o número de acidentes por dia por uma distribuição de Poisson de média lambda = 463.0. Por exemplo, A probabilidade de ter menos do que 400 acidentes por dia é 

> add(evalf(subs(lambda = 463.0, P)), n = 0 .. 400)
0.3170633533e-100 (3.4.1.23)
 

A probabilidade é bem baixa.  

 

Análise de dados 

Regressão polinomial e exponencial 

Se tivermos um conjunto de pontos x[1], y[1], () .. (), x[n], y[n] e queremos achar uma curva com formato pré-determinado que melhor se ajuste a estes pontos, usamos o método dos mínimos quadrados. Por exemplo, suponha que temos o conjunto de pontos 

> L := [[7., 8.3], [2., -.96], [3., -.46], [8., 9.8], [6., 5.5], [10., 12.], [9., 9.8], [1., -2.8], [4., 1.8], [5., .94]]
[[7., 8.3], [2., -.96], [3., -.46], [8., 9.8], [6., 5.5], [10., 12.], [9., 9.8], [1., -2.8], [4., 1.8], [5., .94]] (4.1.1)
 

cujo gráfico é 

> plot(L, style = point, symbolsize = 10)
Plot_2d
 

Observe que os pontos parecem formar uma linha. Para a achar a reta que melhor se ajusta aos pontos, primeiro construimos os seguintes vetores a partir da lista de dados: 

> V1, V2 := Vector(map2(op, 1, L)), Vector(map2(op, 2, L))
Vector[column](%id = 18446744073855372934), Vector[column](%id = 18446744073855373054) (4.1.2)
 

depois usamos o comando PolynomialFit , onde o primeiro argumento é o grau do polinômio: 

> x := 'x'; -1

> reta := PolynomialFit(1, V1, V2, x)
(4.1.3)
 

Vamos sobrepor os gráficos da reta e dos dados para comparação. 

> display({plot(reta, x = 1 .. 10), plot(L, style = point, symbolsize = 10)})
Plot_2d
 

 

Vamos agora analizar um outro conjunto de dados 

> L := [[1., 4.7826], [2., 11.736], [3., 12.986], [4., 25.317], [5., 47.503], [6., 60.753], [7., 115.88], [8., 164.73], [9., 261.67], [10., 469.54]]
L := [[1., 4.7826], [2., 11.736], [3., 12.986], [4., 25.317], [5., 47.503], [6., 60.753], [7., 115.88], [8., 164.73], [9., 261.67], [10., 469.54]]
[[1., 4.7826], [2., 11.736], [3., 12.986], [4., 25.317], [5., 47.503], [6., 60.753], [7., 115.88], [8., 164.73], [9., 261.67], [10., 469.54]]
[[1., 4.7826], [2., 11.736], [3., 12.986], [4., 25.317], [5., 47.503], [6., 60.753], [7., 115.88], [8., 164.73], [9., 261.67], [10., 469.54]]
(4.1.4)
 

Vamos fazer um gráfico dos pontos para verificar alguma tendência. 

> plot(L, style = point, symbolsize = 10)
Plot_2d
 

Uma análise dos pontos mostra que os dados podem ser aproximando por uma função exponencial. 

> V1, V2 := Vector(map2(op, 1, L)), Vector(map2(op, 2, L))
Vector[column](%id = 18446744073878819174), Vector[column](%id = 18446744073878819294) (4.1.5)
 

> func1 := ExponentialFit(V1, V2, x)
`+`(`*`(HFloat(3.499147975656213), `*`(exp(`+`(`*`(HFloat(0.488020149721212), `*`(x))))))) (4.1.6)
 

Caso haja dúvida, por exemplo, alguém pode achar que uma parábola poderia se ajustar melhor, podemos tentar outras regressões. 

> func2 := PolynomialFit(2, V1, V2, x)
`+`(HFloat(81.49068999999989), `-`(`*`(HFloat(55.98848803030303), `*`(x))), `*`(HFloat(8.933396212121213), `*`(`^`(x, 2)))) (4.1.7)
 

Vamos agora fazer o gráfico conjunto. 

> display([plot(L, style = point, symbolsize = 10), plot(func1, x = 1 .. 10, color = red), plot(func2, x = 1 .. 10, color = blue)])
display([plot(L, style = point, symbolsize = 10), plot(func1, x = 1 .. 10, color = red), plot(func2, x = 1 .. 10, color = blue)])
Plot_2d
 

Note que a função de segundo grau tem um péssimo ajuste. 

Método dos mínimos quadrados usando o pacote CurveFitting 

Se tivermos um conjunto de pontos x[1], y[1], () .. (), x[n], y[n] e queremos achar uma curva com formato pré-determinado que melhor se ajuste a estes pontos, usamos o método dos mínimos quadrados. Por exemplo, suponha que temos o conjunto de pontos 

> restart

> with(plots); -1

> L := [[1., -2.8], [2., -.96], [3., -.46], [4., 1.8], [5., .94], [6., 5.5], [7., 8.3], [8., 9.8], [9., 9.8], [10., 12.]]
[[1., -2.8], [2., -.96], [3., -.46], [4., 1.8], [5., .94], [6., 5.5], [7., 8.3], [8., 9.8], [9., 9.8], [10., 12.]] (4.2.1)
 

cujo gráfico é 

> plot(L, style = point, symbolsize = 10)
Plot_2d
 

Queremos descobrir qual é a reta que melhor aproxima estes pontos. O comando LeastSquares do pacote CurveFitting encontra os parâmetros de uma curva que se ajusta ao pontos. Vamos primeiro carregar o pacote. 

A reta é encontrada da seguint maneira: 

> with(CurveFitting); -1

> reta := evalf(LeastSquares(L, x, curve = `+`(`*`(a, `*`(x)), b)), 2)
`+`(`-`(5.1), `*`(1.7, `*`(x))) (4.2.2)
 

O primeiro argumento é a lista de pontos, o segundo é o nome da variável e o terceiro é o formato da curva. No caso acima, usamo a equação da reta. Os gráficos podem ser sobrepostos para melhor vizualização. 

> display({plot(reta, x = 1 .. 10), plot(L, style = point, symbolsize = 10), plot([[7., subs(x = 7., reta)], [7., 8.3]], color = blue, thickness = 4)})
display({plot(reta, x = 1 .. 10), plot(L, style = point, symbolsize = 10), plot([[7., subs(x = 7., reta)], [7., 8.3]], color = blue, thickness = 4)})
Plot_2d
 

Alguns pontos estão mais próximos da reta enquanto outros estão mais afastados. Por exemplo, a distância vertical (veja gráfico) do ponto [7., 8.3] para a reta é 

> `+`(8.3, `-`(subs(x = 7., reta)))
1.5 (4.2.3)
 

Podemos calcular esta distância para cada um dos pontos e elevar cada uma ao quadrado. A reta encontrada pelo método dos mínimos quadrados minimiza a soma das distâncias ao quadrado. No caso acima, o valor mínimo é 

> add(`*`(`^`(`+`(L[i][2], `-`(subs(x = L[i][1], reta))), 2)), i = 1 .. 10)
11.4508 (4.2.4)
 

Qualquer outra reta terá a soma das distâncias ao quadrado maior do que 11.4508. 

 

Exercício 

1(a) Defina uma função f  tal que, dado um valor de x, a função fretorna `+`(`*`(1.7, `*`(x)), `-`(5.4))  mais um valor aleatório com uma distribuição Normal de média 0 e desvio padrão 1. (b) Gere 10 pontos fazendo x  variar de 1 a 10 e faça o gráfico dos pontos. (c) Usando o método dos mínimos quadrados, ache a reta que melhor se ajusta aos pontos do item b. (d) Repita o item c para mais pontos e verifique que a reta obtida se aproxima da reta  

 

2. Faça o ajuste dos pontos [[1., 4.7826], [2., 11.736], [3., 12.986], [4., 25.317], [5., 47.503], [6., 60.753], [7., 115.88], [8., 164.73], [9., 261.67], [10., 469.54]]
[[1., 4.7826], [2., 11.736], [3., 12.986], [4., 25.317], [5., 47.503], [6., 60.753], [7., 115.88], [8., 164.73], [9., 261.67], [10., 469.54]]
[[1., 4.7826], [2., 11.736], [3., 12.986], [4., 25.317], [5., 47.503], [6., 60.753], [7., 115.88], [8., 164.73], [9., 261.67], [10., 469.54]]
[[1., 4.7826], [2., 11.736], [3., 12.986], [4., 25.317], [5., 47.503], [6., 60.753], [7., 115.88], [8., 164.73], [9., 261.67], [10., 469.54]].