# -*- coding: utf-8 -*-
"""
Created on Wed Mar 26 08:41:28 2025

@author: magnu
"""


#Oppgave 1-->4 er lik som Øving 13, så nedenfor starter oppgave 5 med innlimte verdier fra øving 13
    
import numpy as np

#Fastmerker (meter)

HA= 77.237 
HB= 79.582 

#Ukjente er H1, H2

#Observerte høydeforskjeller Δh (meter)
 
ΔhA1=1.916 
Δh12=-0.174 
Δh2B=-0.395 
ΔhB1=0.560 
ΔhA2=1.740 

h_h = np.array([
    ΔhA1,
    Δh12,
    Δh2B,
    ΔhB1,
    ΔhA2
     ])

#Observert avstand D (km)

D_ΔhA1=0.82
D_Δh12=1.22
D_Δh2B=2.03
D_ΔhB1=0.80
D_ΔhA2=1.29

d_h = np.array([
    D_ΔhA1,
    D_Δh12,
    D_Δh2B,
    D_ΔhB1,
    D_ΔhA2
     ])


#Standardavvik til høydeobservasjoner pr km

sD=0.002 + (0.002*d_h)


Std_mean=(sD[0]+sD[1]+sD[2]+sD[3]+sD[4])/5

#Velger apriori

apriori=0.004

#Finne P-matrise

p=apriori**2/sD**2

P=np.array([
    [p[0],0,0,0,0],
    [0,p[1],0,0,0],
    [0,0,p[2],0,0],
    [0,0,0,p[3],0],
    [0,0,0,0,p[4]]
    ])

print("\nFRI UTJEVNING, MANUELL METODE")

A=np.array([
    [1,0,0],
    [-1,1,0],
    [0,-1,1],
    [1,0,-1],
    [0,1,0]
    ])

#Nye frihetsgrader

n=5

e=3

f_fri=n-e


#Pga nye observasjonsligniger blir også l ny

l=np.array([
    [ΔhA1+HA],
    [Δh12],
    [Δh2B],
    [ΔhB1],
    [ΔhA2+HA]
    ])


#Finne ny Normalligningsmatrisen

N=A.T@P@A

#Finne Normalligningens høyre side

h=A.T@P@l


#Estimere høydene til H1 og H2 med minste kvadraters metode

x=np.linalg.inv(N)@h

ny_H1_hatt=x[0]
ny_H2_hatt=x[1]
HB_hatt=x[2]


#Finne residualer, v

v=A@x-l

#Feilkvadratsum, V^TPV

Feilkvadratsum=v.T@P@v

#Stadardavviket til vektsenheten p=1

S0=np.sqrt(Feilkvadratsum/f_fri)

#Foreta en observasjonstest basert på multippel t-testing.
#i. For hver observasjon; estimer sannsynlig grovfeil, ∇, med tilhørende
#standardavvik, 𝑠∇, og testverdi, T.

#For å teste observasjonene må vi sette opp ny designmatrise A_test

A_test=np.array([
    [1,0,0,0],
    [-1,1,0,0],
    [0,-1,1,0],
    [1,0,-1,0],
    [0,1,0,1]
])

#Ny N_test for testing av observasjoner
N_test=A_test.T@P@A_test


#Finne Normalligningens nye høyre side, h_test1

h_test=A_test.T@P@l

#Bruker minste kvadraters metode for teste første observasjon
    
x_test=np.linalg.inv(N_test)@h_test

print(f"\nGrovfeil, ∇ er:{x_test[3]})")

#Beregner standardavviket til grovfeil, ∇

#Finne nye residualer, v

v_test=A_test@x_test-l

Feilkvadratsum_test=v_test.T@P@v_test


#Regner ut ny S0
    #Èn mindre frihetsgrad nå som vi har en ukjent ekstra

S0_test=np.sqrt(Feilkvadratsum_test/(f_fri-1))


#Finne ny Qxx_test
    
Qxx_test=np.linalg.inv(N_test)


#Standardavviket til grovfeilen til den første observasjonen

sD_gf=S0_test*np.sqrt(Qxx_test[3,3])

print(f"\nStandardavviket til grovfeilen, ∇ er:{sD_gf}")

#Finne T-verdien til første grovfeil

T=x_test[3]/sD_gf

print(f"\nT-verdien til grovfeilen er: {T}")\


    #Multippel T-test

#Observasjon 1
    # Grovfeil  ∇1 er: 0.00455662
    # Std.avvik ∇1 er: 0.00751144
    # T-verdien ∇1 er: 0.60662463

#Observasjon 2
    # Grovfeil  ∇2 er: -0.00248675
    # Std.avvik ∇2 er: 0.00759283
    # T-verdien ∇2 er: 0.32752034

#Observasjon 3
    # Grovfeil  ∇3 er: -0.0097313
    # Std.avvik ∇3 er: 0.00214756
    # T-verdien ∇3 er: 4.53117273

#Observasjon 4
    # Grovfeil  ∇4 er: -0.00973095
    # Std.avvik ∇4 er: 0.00214756
    # T-verdien ∇4 er: 4.53117273

#Observasjon 5
    # Grovfeil  ∇5 er: -0.00455662
    # Std.avvik ∇5 er: 0.00751144
    # T-verdien ∇5 er: 0.60662463
    

#Tester om mest sannsynlige grovfeil er signifikant forskjellig fra 0

#Test av signifikansnivå for enkelttester:

alfaTot = 0.05 #ønsket signifikansnivå

alfa_im = 1 - (1-alfaTot)**(1/n)    #_im står for alfa i "MERKET"
alfa_i = alfa_im / 2            #justerer for to-sidig test
print(f"\nAlfa_i er:{alfa_i}")
print("\nSignifikansnivå for observasjon i:") #Sammenlikner man T-verdiene med denne verdien eller
print(alfa_i)               #bruker man denne verdien til å slå opp i t-tabellen?
                            #Hvilken er den kritiske verdien? alfa_i eller t_tabell?
from scipy.stats import t

t_tabell = t.ppf(1 - alfa_i, n-e-1)
print("\nfra t-tabell:")
print(t_tabell)


    #Gjennomfører grunnlagstest!

#Indre pålitelighet

I=np.array([            #Identitetsmatrise (nxn-dimensjon)
    [1,0,0,0,0],
    [0,1,0,0,0],
    [0,0,1,0,0],
    [0,0,0,1,0],
    [0,0,0,0,1]
     ])

Qll=np.linalg.inv(P)

F=A@np.linalg.inv(N)@A.T@P-I

v_R=F*l

Qvv=F@Qll@F.T

R=Qvv@P

print(f"\nRedundansmatrisen R er:\n{R}")

    #Limer inn v^TPv fra tvangsutjevning i oppgave 13!
    
Feil_tvang=1.09866893

delta_f=1

F_test=((Feil_tvang-Feilkvadratsum)/delta_f)/((Feilkvadratsum/f_fri))

f_tvang=3

from scipy.stats import f

# Frihetsgrader
df1 = f_tvang - f_fri # Teller
df2 = f_fri     # Nevner

# Kritisk F-verdi for alpha = 0.05
alpha = 0.05
F_kritisk = f.ppf(1 - alpha, df1, df2)

# Sammenligning
print(f"F-observert: {F_test}")
print(f"F-kritisk (5% signifikansnivå): {F_kritisk:.4f}")

if F_test > F_kritisk:
    print("Det ER signifikant tvang i grunnlaget!")
else:
    print("Det er IKKE signifikant tvang i grunnlaget")