Vi introducerar en familj av implicita probabilistiska integratörer för initiala värdeproblem (IVPs), med utgångspunkt från flerstegsmetoden Adams-Moulton. Den implicita konstruktionen möjliggör dynamisk återkoppling från det kommande tidssteget, i motsats till tidigare probabilistiska integratörer, som alla bygger på explicita metoder. Vi börjar med en kortfattad undersökning av det snabbt expanderande området för probabilistiska ODE-lösare. Vi introducerar sedan vår metod, som bygger på och anpassar Conrad et als arbete. (2016) och Teymur et al. (2016), och ge ett rigoröst bevis på dess väldefinieradhet och konvergens. Vi diskuterar problemet med kalibreringen av sådana integratörer och föreslår ett tillvägagångssätt. Vi ger ett belysande exempel som belyser effekten av användningen av probabilistiska integratörer-inklusive vår nya metod-i inställningen av parameter inferens inom ett omvänt problem. Vi anser att det gemensamma statistiska problemet med att härleda modellparametrar från data Y. I en Bayesian inställning, parametern bakre ges av p( skaY )  på p(Y  på)p( ska). Antag att vi har en regressionsmodell där sannolikhetsbegreppet p (Y) kräver att vi löser en vanlig differentialekvation (ODE). För varje datum har vi Y j = x(t Yj ) + ε j för någon latent funktion x(t) som uppfyller = f (x, ) och vektor för mätfel ε med spridningsparameter σ. Vi kan skriva den fullständiga modellen som p( Eftersom x är latent, ingår det som en integrerad del av den bakre modellen. Denna mer allmänna nedbrytning skulle inte behöva beaktas explicit i t.ex. en linjär regressionsmodell för vilken x =  på 1 + 2 t; här skulle vi helt enkelt ha p(x på) = δ x ( på 1 + på 2 t). Med andra ord, med tanke på att det inte finns någon osäkerhet i x, och modellen skulle minska till helt enkelt p( ska, σ skaY ) p(Y  på, σ)p( ska)p(σ). I vårt fall, dock, x definieras implicit genom ODE till = f (x,  på) och p (x på) är därför inte längre trivialt. Vad vi menar med detta är att x bara kan beräknas ungefär och därmed följa den centrala principen för probabilistiska numeriska metoder (Hennig et al., 2015)-vi tilldelar det en sannolikhetsfördelning som representerar vår brist på kunskap om dess sanna värde. Vårt fokus här är på initiala värdeproblem (IVPs) där vi antar att det initiala värdet X 0 x(0) är känt (även om en förlängning till okänt X 0 är okomplicerat). Vi har alltså installationen p( ska, σ, x skay, X 0 ) p(Y  skax, σ)p(x ska, X 0 )p( ska)p(σ). (1) För våra syften, den intressanta termen på höger sida är p(x på), där vi nedan utelämnar X 0. I vid bemärkelse är vårt mål att så exakt som möjligt redogöra för det numeriska fel som är oundvikligt i beräkningen av x, och att göra detta inom en probabilistisk ram genom att beskriva p(x på). Vi vill sedan ta hänsyn till effekten av denna osäkerhet, eftersom den sprids genom (1), när man gör en slutledning i en omvänt problematik. Ett experimentellt exempel i detta sammanhang behandlas i avsnitt 3.  Innan vi ger en sammanfattning av det nuvarande tillståndet av probabilistiska numeriska metoder (PN) för ODEs, tar vi en kort avledning. Det är intressant att notera att begreppet att definiera en distribution för p(x) har framträtt i den senaste litteraturen i olika former. Till exempel, en serie papper (Calderhead et al., 2009; Dondelinger et al., 2013; Wang och Barber, 2014; Macdonald et al., 2015), som uppstod separat från PN i sin moderna form, försöker undvika att lösa ODE helt och hållet och istället ersätta den med en "surrogat" statistisk modell, parameteriserad av  och med den primära motivationen att minska den totala beräkningen. Den centrala principen i dessa papper är att utföra gradientmatchning (GM) mellan de fullständiga och surrogatmodellerna. En följd av denna ram är införandet av en distribution p (x på,  och) liknande p (x på) som visas i (1). Syftet är att uppskatta x med hjälp av statistiska tekniker, men det finns inget försök att modellera felet i sig, helt enkelt ett försök att minimera skillnaden mellan den sanna lösningen och dess surrogat. Dessutom är parametrarna i de surrogatmodeller som föreslås i GM-ramverket monterade genom konditionering på data Y, vilket betyder att p(x på, ) måste betraktas som en datakonditionerad bakre p(x på,  och Y ). Enligt vår mening är detta problematiskt, eftersom när osäkerheten i en kvantitet av intresse uppstår enbart på grund av att de numeriska metoder som används för att beräkna den är inexakta, bör slutsatsen om denna kvantitet inte baseras på uppgifter som är resultatet av ett experiment. Den cirkuläritet som induceras i (1) av Y-konditionering är klar. Den grundläggande förändringen i tänkandet i tidningarna av Hennig och Hauberg (2014) och Chkrebtii et al. (2016), bygger på Skilling (1991), och sedan följt upp och utvidgats av Schober et al. (2014), Conrad m.fl. (2016), Teymur et al. (2016), Kersting och Hennig (2016), Schober et al. (2018) och andra är det av vad som utgör "data" i den algoritm som används för att bestämma p(x). Till skillnad från GM-metoden används inte de experimentella uppgifterna Y för att konstruera denna distribution. Även om det redan har sagts underförstått i vissa av dessa verk, hävdar vi att detta utgör den viktigaste skillnaden i filosofi. I stället bör vi sträva efter att kvantifiera den numeriska osäkerheten i x först, och sedan sprida denna osäkerhet via datasannolikheten till Bayesian inversion som används för att dra slutsatsen. Detta är i praktiken direkt probabilistisk modellering av det numeriska felet och är den metod som används i PN. Hur kan då x härledas i PN? Den gemensamma tråden här är att en diskret väg Z  på Z 1:N genereras som numeriskt approximerar X  på X 1:N -den diskretiserade versionen av den sanna lösningen x-då "modell förhör" (Chkrebtii et al., 2016 ) F:= f (Z,  på) är tänkt som en typ av numeriska data och x härleds utifrån dessa. Gjort på detta sätt, en helt modellbaserad beskrivning av osäkerheten i x resultat, utan användning av experimentella data Y. En annan central egenskap hos PN-lösare från Chkrebtii et al. (2016) Vidare är att behandla problemet sekventiellt, på ett sätt som en klassisk IVP integrator. I alla GM papper, och faktiskt i Hennig och Hauberg (2014), X behandlas som ett block - härleds alla på en gång, med tanke på data Y (eller, i Hennig och Hauberg, F ). Detta begränsar nödvändigtvis graden av återkoppling som är möjlig från dynamiken i den faktiska ODE, och i en allmän IVP kan detta vara källan till betydande inkorrekthet, eftersom fel i de inexakta värdena Z approximerande X förstärks av själva ODE. I ett sekventiellt tillvägagångssätt, är de numeriska data inte ett statiskt redan existerande objekt som de verkliga data Y är, utan snarare genereras när vi går genom att upprepade gånger utvärdera ODE på en sekvens av indata koordinater. Således är det tydligt att de numeriska data som genereras vid tidpunkten t påverkas av den infererade lösningen ibland före t. Denna iterativ information återkoppling är kvalitativt mycket mer som en standard IVP lösare än en block inference metod och liknar principen för statistisk filtrering (Särkkä, 2013). Vi granskar nu de befintliga dokumenten på detta område närmare för att ge ett sammanhang till vårt eget bidrag i nästa avsnitt. I Chkrebtii et al. (2016) en Gaussian process (GP) tidigare placeras gemensamt över x och dess derivat, sedan i steg i de nuvarande GP parametrar används för att förutsäga ett värde för staten i nästa steg, Z i+1. Detta omvandlas sedan för att ge Modelleringsantagande som nu gjorts är att detta värde fördelas runt det verkliga värdet av derivatan på i+1 med Gaussian fel. På grundval av detta assimileras det nya datumet i modellen, vilket ger en bakre för (x, på) som kan användas som start före i nästa steg. Detta tillvägagångssätt gör inte direkt användning av sekvensen Z; snarare är det bara genereras för att producera de numeriska data F som sedan jämförs med den tidigare modellen i derivatat utrymme. Resultatet är en fördelnings Gaussian bakre över x överensstämmer med sekvensen F. 2 Conrad et al. (2016) inta en annan hållning. Bearbetning av problemet i en diskret inställning, de producerar en sekvens Z av värden som approximerar X, med F i+1 f (Z i+1,  på) utgör data och Z i+1 beräknas från de tidigare värdena Z i och F i genom att använda en viss iterativ relation liknar en randomiserad version av en standard IVP lösare. Observera att det inte finns något försök att kontinuerligt assimilera de genererade värdena till modellen för det okända X eller till och med under körningen av algoritmen. I stället är motiveringen för metoden post hoc i form av ett konvergensteorem som begränsar den maximala förväntade kvadrat-terror max i E på X i 2. En förlängning till flerstegsmetoder-i vilka Z i+1 är tillåtet att bero på flera tidigare värden F ≤i -införs i Teymur et al. (2016) och använder samma grundläggande metod. Olika utökningar och generaliseringar av de teoretiska resultaten i dessa dokument ges i Lie et al. (2017), och en relaterad idé där stegstorleken är randomiserad föreslås av Abdulle och Garegnani (2018). Detta tillvägagångssätt är intuitivt, vilket möjliggör modifierade versioner av standardalgoritmer som ärver kända användbara egenskaper, och ger bevisbara förväntade felgränser. Det är också mer allmänt eftersom det tillåter icke-parametriska bakre fördelningar för x, även om det bygger på Monte Carlo provtagning för att ge empiriska approximationer till det. Mathematically, we write Here, Z  på Z 1:N är approximationen till den okända diskretised lösning funktion X X 1:N, och varje Använda terminologin i Hennig et al. (2015) motsvarar de två beståndsdelarna i den teleskopiska nedbrytningen på höger sida av (2) "beslutsregeln" (hur algoritmen genererar en ny datapunkt F i ) respektive "generativ modell" (som kodar sannolikhetsmodellen för Z). Observera att ur statistisk synvinkel definierar metoden uttryckligen en fördelning över numeriska lösningar Z snarare än en osäkerhet som är centrerad kring den sanna lösningen x (eller X). Konvergensanalysen garanterar sedan förhållandet mellan åtgärden över Z och åtgärden över X. Termen p(F i  påZ i,  på) är tagen i både Conrad et al. (2016) och Teymur et al. (2016) Att vara helt enkelt en deterministisk omvandling; detta kan skrivas i distributionsform som ges av Conrad et al. som en Gaussian centrerad kring utdata Z det i+1 av någon deterministisk singelsteg IVP lösare, med varians skalad i enlighet med begränsningarna i deras teorem. Teymur m.fl. införa en konstruktion för denna term som tillåter konditionering på flera tidigare F i s och som har ett medelvärde som motsvarar flerstegsmetoden Adams-Bashforth. De ger motsvarande allmänna konvergensresultat. Deras bevis är också lätt verifieras för att vara giltig för implicita flerstegsmetoder-ett resultat vi vädjar till senare-även om den specifika implicita integrator modell de föreslår är metodologiskt inkonsekvent, av skäl som vi kommer att förklara i Avsnitt 2. I alla de metoder som beskrivits hittills, Monte Carlo provtagning krävs för att marginalisera F och därmed beräkna p (Z på). Detta utgör en märkbar beräkningsgrund. Ett tredje tillvägagångssätt, relaterat till stokastisk filtrering, presenteras i Schober et al. (2014), Kersting och Hennig (2016) och Schober et al. (2018). Dessa dokument utvecklar ett ramverk som inte bygger på provtagning, utan istället gör det enklare antagandet att alla distributioner är Gaussiska, och sprider osäkerheten från steg till steg med hjälp av teorin om Kalmanfiltrering (Särkkä, 2013). Detta är en alternativ lösning till avvägningen mellan exakthet och beräkningskostnad, en punkt som erkänns i dessa dokument. För jämförelsens skull kan vi löst skriva om deras allmänna tillvägagångssätt i vår notation enligt följande: där vi skriver x [i] istället för Z i för att betona att detta representerar en i-times uppdaterad modell för den kontinuerliga lösningen x, snarare än i'th iteration av en algoritm som genererar en approximation till diskret X i. Denna algoritm förutspår ett värde för den nya stateZ i+1 från den nuvarande modellen och alla tidigare data, sedan genererar en datapunkt baserad på den förutsägelsen, och uppdaterar sedan modellen baserat på detta nya datum. Observera att alla distributioner inom detta ramverk är Gaussiska, för att tillåta snabb filtrering, och som ett resultat icke-linjäriteterna i f är effektivt linjäriserade, och varje numerisk metod som producerar icke-gaussiska fel har Gaussians monteras på dem ändå. Denna filtreringsmetod är intressant på grund av det tidigare angivna nedfallet av maximering av graden av återkoppling från ODE-dynamiken till lösaren. Den prognostiserade utvärderingsuppdateringen 3 som föreslås av (3) innebär att information från ODE-funktionen vid nästa steg t i+1 matas tillbaka till förfarandet vid varje steg, till skillnad från andra metoder som bara förutsäger framåt. I numerisk analys är detta typiskt en efterlängtad egenskap, vilket leder till metoder med förbättrad stabilitet och noggrannhet. Det är dock fortfarande en tredelad procedur, liknande till exempel att para AdamsBashforth och Adams-Moulton integratörer som används i PEC-läge (Butcher, 2008). Detta samband hänvisas till i Schober et al. (2018). Vi föreslår nu ett annat, nytt sekventiellt förfarande som också innehåller information från ODE vid tidpunkten för steg t i+1 men gör det direkt. Detta producerar en sann implicit probabilistisk integratör, utan en subtil inkonsekvens som finns i den metod som föreslås av Teymur et al. (2016). Där definierar analogen (2) en gemensam Gaussisk fördelning över Z i+1 och F i+1 (den högra komponenten, med F ≤i ersatt av F ≤(i+1) ) men genererar sedan F i+1 genom att passera Z i+1 genom funktionen f (den vänstra komponenten). Detta ger två ömsesidigt inkompatibla betydelser till F i+1, en linjärt och en icke-linjärt relaterad till Z i+1. Vår föreslagna metod löser detta problem. Vi utnyttjar just skillnaden i dessa två kvantiteter genom att skilja dem åt och direkt bestraffa skillnaden mellan dem. För att introducera idén betraktar vi först det endimensionella fallet, sedan generaliserar vi till ett flerdimensionellt sammanhang. Vi noterar först att till skillnad från i de explicita randomiserade integratörerna av Conrad et al. (2016) och Teymur et al. (2016) Vi har inte tillgång till den exakta deterministiska Adams-Moulton prediktorn, som vi sedan skulle kunna lägga till en noll-metal perturbation. Det krävs därför ett alternativt tillvägagångssätt. Betrakta istället följande fördelning som direkt avancerar integratorn ett steg och beror endast den aktuella punkten: Här, r(z) är en positiv diskrepansmått i derivatutrymme definieras i nästa stycke, och g är en η-skalad funktionell omvandling som säkerställer att uttrycket är en giltig sannolikhetsfördelning i variabeln z. Ett konkret exempel belyser definitionen. Betrakta den enklaste implicita metoden, bakåt Euler. Detta definieras av förhållandet Z i+1 = Z i + hF i+1 och kan normalt endast lösas genom en iterativ beräkning, eftersom F i+1  på (Z i+1,  på) naturligtvis är okänd. Om den slumpmässiga variabeln Z i+1 har värdet z, då kan vi uttrycka F i+1 som en funktion av z. Specifikt har vi F i+1 (z) = h −1 (z − Z i ). Skillnaden r(z) mellan värdet av F i+1 (z) och värdet av f (z,  på) kan sedan användas som ett mått på felet i den linjära metoden, och bestraffas. Detta är detsamma som att bestraffa skillnaden mellan de två olika uttrycken för F i+1 som härrör från den tidigare beskrivna naiva utvidgningen av (2) till det implicita fallet. Vi skriver Comparing (4) och (5), r(z) är uttrycket h, och g är i detta fall omvandling u → exp(−u 2 /2η 2 ). Detta tillvägagångssätt främjar direkt lösaren i ett enda steg, utan att samla in explicita numeriska data som i tidigare tillvägagångssätt. Den är i allmänhet icke-parametrisk och kräver antingen provtagning eller approximation för att vara användbar (mer om vilken i nästa avsnitt). Eftersom f i allmänhet är icke-linjär, följer det att r är icke-linjär också. Därav följer att densiteten i ekvation (5) inte resulterar i ett Gaussiskt mått, trots att g är en kvadrat-exponentiell omvandling. Generaliseringen till högre order implicita linjära flerstegsmetoder av Adams-Moulton (AM) typ, med formen Följande analys visar den väldefinierade och konvergens egenskaper av konstruktionen som föreslås i avsnitt 2. Först visar vi att fördelningen (6) är väldefinierad och korrekt, genom att bevisa 4 ändlighet och strikt positivitet av den normaliserande konstanten K. Vi beskriver sedan villkoren på h-beroende av skalningsparametern η, så att de slumpmässiga variablerna  och i (8) uppfyller hypoteserna i Teorem 3 i Teymur et al. (2016). Konvergensen i vår metod beror särskilt på Adams-Moulton-analogen av detta resultat. Redovisas av  till h,s : R d×s → R d den deterministiska kartan definieras av s-steg Adams-Moulton metoden. Till exempel är den implicita kartan som hör ihop med den bakåtvända Euler-metoden – "nollstegsmetoden" AM-metoden – för en fast parameter  ska vara.h  och,0 (Z i ) = Z i + hf (- h -,0 (Z i ), -). Mer allmänt, den karta som förknippas med s-steg AM metoden är där, och β j  på R + är Adams-Moulton koefficienter. Observera att den deterministiska Adams-Moulton uppskattningen för Z i+1 (Z i‐s+1:i) representerar den deterministiska Adams-Moulton uppskattningen för Z i+1. Med tanke på ett sannolikhetsutrymme (., F, P), definiera för varje i  på N den slumpmässiga variabeln  för förhållandet mellan uttrycken (6) och (8) behandlas i del (i) av följande Theorem, vars bevis ges i det kompletterande material som åtföljer detta papper. Teorem. Antag att vektorfältet f (·,  på) är globalt Lipschitz med Lipschitz konstant L f, på > 0. Rätta s  och R q, och 0 < h < (L f, och β −1 ) −1. Om η = kh ρ för vissa k > 0 oberoende av h och ρ ≥ −1, då följande uttalanden innehar: i) Den funktion som definieras i (6) är en väldefinierad sannolikhetstäthet. (ii) För varje r ≥ 1, det finns en konstant 0 < C r < ∞ som inte beror på h, sådan att, Den probabilistiska integratorn definieras av (6) konvergerar i medel-kvadrat som h → 0, i samma takt som den deterministiska s-steg Adams-Moulton metoden. Extrapoleringsdelen av en linjär metod fungerar separat på varje komponent i ett flerdimensionellt problem. i−j för varje komponent k i tur och ordning. Naturligtvis är detta inte sant om omvandlingen Z i+1 → F i+1 f (Z i+1 ), utom i det triviala fallet där f är linjär i z; alltså i (2) den högra fördelningen är komponentvis oberoende medan den vänstra är inte. Alla tidigare sekventiella PN-integratörer har behandlat det flerdimensionella problemet på detta sätt, som en produkt av endimensionella relationer. I vårt förslag är det inte meningsfullt att överväga systemet med ekvationer komponent för komponent, på grund av förekomsten av den icke-linjära f (z,  på) term, som visas som en inneboende del av den steg framåt fördelning p(Z i+1  och Z ≤i,  och η). Den multidimensionella analogen (6) bör ta hänsyn till detta och definieras över alla d-mått tillsammans. För vektorvärderad z, Z k, F k, definierar vi därför var r(z) = β är nu en d × 1 vektor av avvikelser i derivatutrymme, och H är en d × d matris som kodar lösaren skala, generaliserande η. Rättframma ändringar av beviset ger multidimensionella analoger till påståendena i Theorem. Frågan om kalibrering av ODE-lösare behandlas utan konsensus i varje behandling av detta ämne som nämns i avsnitt 1. Tillvägagångssätten kan i stort delas upp i de av "framåt"-typ, där det finns ett försök att direkt modellera vad den teoretiska osäkerheten i ett lösarsteg bör vara och propagera det genom beräkningen; och de av "bakåt"-typ, där osäkerhetsskalan på något sätt matchas efter beräkningen till den som föreslås av någon annan indikator. Båda dessa brister, de förra på grund av den inneboende svårigheten att uttryckligen beskriva felet, och de senare på grund av att det per definition är mindre exakt. En viktig stötesten är att det i allmänhet är ett utmanande problem att till och med definiera vad det innebär för en uppskattning av osäkerheten att vara väl kalibrerad. I detta dokument kräver vi ett sätt att sätta H. Vi går vidare genom att modifiera och generalisera en idé från Conrad et al. (2016) som hör till kategorin "bakåt". Där, variansen av steg framåt fördelning Var(Z i+1  på · · · ) anses vara en matris Z = αh ρ I d, med α bestäms av en skala-matchning förfarande som säkerställer integratorns utgångar en global felskala i linje med förväntningarna. Vi hänvisar läsaren till den detaljerade beskrivningen av detta förfarande i avsnitt 3.1 i det dokumentet. Konvergensen är dessutom ett resultat av Teymur et al. (2016) innebär att exponenten ρ för det probabilistiska s-steget Adams-Bashforth integrator bör tas till 2s + 1. I vår metod, vi inte kan relatera en sådan matris Z direkt till H eftersom från definitionen (9) är det tydligt att H är en skalmatris för spridningen av derivatan F i+1, medan Z mäter spridningen av staten Z i+1. För att omvandla till rätt utrymme utan att linearisera ODE, tillämpar vi multivariat delta-metoden (Oehlert, 1992) för att ge en approximation för variansen av den transformerade slumpvariabeln, och ställa H att vara lika med resultatet. Således där J f är Jacobian av f. Medelvärdet E(Z i+1 ) är okänt, men vi kan använda en explicit metod av lika eller högre ordning för att beräkna en uppskattning Z AB i+1 till försumbar kostnad, och använda J f (Z AB i+1 ) istället, under antagandet att dessa är rimligt nära. Kom ihåg att vi ungefär kalibrerar metoden så att en viss nivå av approximation är oundviklig. Denna kommentar gäller även i de fall där jakobinen inte är analytiskt tillgänglig och beräknas numeriskt. Sådana approximationer påverkar inte algoritmens grundläggande konvergensegenskaper, eftersom de inte påverkar den stegvisa fördelningens h-skalning. Vi noterar också att vi bara matchar varianser / spridda parametrar och ingenstans förutsatt att fördelningen (9) är Gaussian. Denna idé har en viss likhet med det ursprungliga konceptet i Skilling (1993), där en skalär "stelhetskonstant" används på ett liknande sätt för att omvandla osäkerhetsskalan från lösningsutrymme till derivatutrymme. Vi fastställer nu lämplig h-skalning för H genom att ställa in exponenten ρ. Det villkor som krävs av den unika analysen i detta dokument är att η = kh ρ; del (iii) av Theorem visar att vi behöver ρ ≥ s + för en h-oberoende positiv-definit matris Q. Eftersom J f är oberoende av h, betyder detta att vi måste ställa in Z att vara proportionell mot h 2 s+ 1 2 ), och därmed har vi Vår konstruktion har den positiva följden av att ge en icke-trivial korskorrelation struktur till felkalibreringsmatrisen H, vilket möjliggör en rikare beskrivning av felet i flerdimensionella problem, något frånvarande från tidigare metoder. Dessutom härleds denna ytterligare information via direkt återkoppling från ODE, som vi har visat är ett önskvärt attribut. I den form som beskrivs i föregående avsnitt resulterar vår algoritm i en icke-parametrisk fördelning för Z i+1 vid varje steg. Med detta tillvägagångssätt kan en beskrivning av osäkerheten i den numeriska metoden endast utvärderas genom ett provtagningsförfarande enligt Monte Carlo vid varje iteration. Även om denna provtagning utförs med hjälp av en metod som är väl lämpad att rikta distributioner nära Gaussian-vi använder en modifierad version av den förkonditionerade Crank-Nikolson algoritm föreslås av Cotter et al. (2013) – det finns tydligt en betydande beräkningspåföljd i samband med detta. Det enda sättet att undvika denna påföljd är genom att återgå till distributioner av standardformulär, som är lätta att ta prover från. En möjlighet är att ungefär (6) av en Gaussisk distribution-beroende på hur denna approximation utförs desideratum av att upprätthålla information återkoppling från den framtida dynamiken i målfunktionen kan bibehållas. Till exempel, en första ordning Taylor expansion av, När substitueras till r(z) enligt definitionen i (9), ger en catimatioñ r(z) som är linjär i z. Detta ger en icke-centrerad Gaussian när den omvandlas till en sannolikhetsmått som i (9). Definiera (hβ Vi tar detta tillfälle i akt att påminna läsaren om den olyckliga konventionen från numerisk analys som resulterar i att s har olika betydelser här och i föregående stycke - den explicita metoden av ordning s är den med s steg, medan den implicita metoden av ordning s är den med s - 1 steg. 6 några enkla algebra ger ögonblicken av den approximerande Gaussian åtgärd för nästa steg som μ = Z i + ρ −1 w och Vi noterar att detta förfarande är bara för att underlätta enkel provtagning-tanker(z) är linjär i z, införandet av den första ytterligare termen från Taylor expansion innebär att information om icke-linjäritet (i z) f fortfarande införlivas till andra ordningen, och den genererade lösningen Z är inte gemensamt Gaussian över tid steg i. Eftersom punkt 1 är ordning 1 i h påverkar denna approximation inte den globala konvergensen för integratorn, så länge som H fastställs i enlighet med de principer som beskrivs i punkt 2.3. Denna metod för att lösa implicita integratorer genom att linjärisera dem i f är välkänd i klassisk numerisk analys, och de resulterande metoderna kallas ibland semiimplicit metoder (Press et al., 2007). Vi illustrerar vår nya algoritm med tanke på fallet med ett enkelt omvänt problem, FitzHughNagumo modellen diskuteras i Ramsay et al. (2007) och därefter behandlas i ett antal dokument i denna fråga. Detta är ett tvådimensionellt, icke-linjärt dynamiskt system med tre parametrar  och = (en, två, tre, tre), vars värden (en = 0,2, två = 0,2, tre = 3,0) väljs för att producera periodisk rörelse. Med problemet att ha en Bayesian struktur, skriver vi ner den bakre som Detta uttryck erinrar (1), men med Z ersätta x som beskrivs i avsnitt 1.2. För enkelhetens skull använder vi det kända värdet av σ i hela, så inkludera det inte i den bakre modellen. Conrad och Al. (2016) anmärkning om den fördomsfulla effekten på den bakre fördelningen av den naiva utvärderingen av framåtmodellen med hjälp av en standardnumerisk metod. De visade att deras probabilistiska integratör återgäldar större posteriors, vilket förhindrar missriktat övertro i en felaktig uppskattning. Vi utvidgar nu dessa experiment till vår nya metod. I förbigående noterar vi intressant aktuell teoretisk utveckling som diskuterar den kvantitativa effekten på den bakre slutsatsen av randomiserade framåtmodeller, som presenteras i Lie et al. (2018).  Den första uppgiften är att härleda det lämpliga värdet för den totala skalningskonstanten α för varje metod, som ska användas vid inställning av matrisen H i (10). Som i Conrad et al. (2016), vi beräknar ett värde α * som maximerar överenskommelsen mellan utfallet av den probabilistiska integratorn och ett mått på globala fel från den deterministiska metoden, och sedan fixa och gå vidare med detta värde. För var och en av flera metoder M, α * M beräknades för ett intervall av värden av h och var nära konstant genom hela, vilket tyder på att h-skalning förespråkas i avsnitt 2.3 (dvs.. ta jämlikhet i den bundna delen (iii) av Theorem) är den korrekta. Denna punkt har inte specifikt tagits upp i tidigare arbeten i denna fråga. Den faktiska maxima α * M för varje metod är olika och ytterligare forskning krävs för att undersöka om ett samband kan härledas mellan dessa värden och någon känd egenskap hos varje metod, såsom antal steg s eller lokal fel konstant av den underliggande metoden. Dessutom förväntar vi oss att dessa värden ska vara problemberoende. I detta fall fann vi α * Efter att ha kalibrerat den probabilistiska integratorn, illustrerar vi dess typiska utgång i figur 1: den översta rutan ritar vägen för 500 iterationer av den probabilistiska bakåt Euler metoden kör vid α = α * AM0 = 0,2. Vi ritar de diskreta värdena Z 1:N för varje repetition, utan att försöka skilja banorna från enskilda körningar. Detta är att betona att varje randomiserad körning (till följd av en annan instantiering av  på) inte är avsedd att ses som ett "typiskt prov" från vissa underliggande kontinuerliga sannolikhetsmått, som i vissa andra probabilistiska ODE-metoder, utan snarare att de tillsammans bildar en ensemble från vilken en empirisk fördelning som karakteriserar diskret-tid lösning osäkerhet kan beräknas. Den nedre rutan ritar samma data men med skuggade band som representerar intervallerna 1σ, 2σ och 3σ, och en streckad linje som representerar det empiriska medelvärdet. Vi överväger nu det omvända problemet med att härleda parametrarna för Fitzhugh-Nagumo-modellen i intervallet t  och [0, 20]. Vi genererar först syntetiska data Y ; 20 tvådimensionella data-punkter samlas in vid tider t Y = 1, 2,..., 20 korrumperas av centrerad Gaussian brus med varians σ = (0.01) · I 2. Vi behandlar sedan parametrarna som okända och kör en MCMC algoritm-Adaptive Metropolis Hastings (Haario et al., 2001 )-för att dra slutsatsen deras bakre distribution. I Conrad och Al. (2016), den ekvivalenta algoritmen utför flera upprepningar av framåt lösa vid varje steg av den yttre MCMC (var och en med en annan omedelbarion av  på) marginaliserar sedan ut för att bilda en förväntad sannolikhet. Detta är beräkningsmässigt mycket dyrt; i våra experiment finner vi att för MCMC att blanda väl, många tiotals upprepningar av framåt lösning krävs vid varje steg. I stället använder vi ett Metropolis-inom-Gibbs-system där vid MCMC iteration k föreslås och godkänns och avvisas att ha fått sin sannolikhet beräknad med samma prov 0:N sedan kan provtas och sannolikhetsvärdet omberäknas redo för nästa förslag. Förslaget i steg k + 1 jämförs sedan med detta nya värde. Pseudo-kod för denna algoritm anges i det kompletterande materialet. Vårt tillvägagångssätt kräver att p(Z på,  på) räknas om exakt en gång för varje gång en ny parameter värde  på * accepteras. Kostnaden för denna strategi begränsas därför av två gånger kostnaden för en MCMC-anläggning som drivs med en deterministisk integratör, som endast uppnås i det scenario som alla föreslagna åtgärder är accepterade. Algoritmen, i motsats till kalibreringsproceduren (som är relativt kostsam men som bara behöver utföras en gång), har därför begränsat ytterligare beräkningsmässiga omkostnader jämfört med den naiva metoden med en klassisk metod. Figur 2 visar kärndensitetsuppskattningar som ungefär motsvarar den bakre fördelningen för de framåtgående Euler, probabilistiska framåt Euler, bakåt Euler och probabilistiska bakåtvända Euler-metoderna. Var och en representerar 1000 parameterprover från simuleringar som körs med stegstorlekar h = 0,005, 0,01, 0,02, 0,05. Detta görs av 11000 totala prover, med de första 1000 kasserade som inbrännning, och resten tunnas ut med en faktor på 10. För varje metod M, dess förberäknade kalibreringsparameter α * M används för att ställa in variansen av  till. Vid större stegstorlekar ger deterministiska metoder både överkonfidentiella och partiska uppskattningar (på olika sidor av det verkliga värdet). I enlighet med resultaten från Conrad et al. (2016) Nederst till höger visar samma effekt med den probabilistiska bakåtvända Euler-metoden som vi har infört i detta dokument. Vi hittar liknande resultat för andra och högre ordningsmetoder, både explicita och implicita. Effektens omfattning är dock relativt liten när det gäller ett så enkelt testproblem, där en högre grad av integrator inte förväntas ge upphov till mycket fel i den framtida lösningen. Vidare arbete kommer att undersöka tillämpningen av dessa metoder på mer utmanande problem. I detta dokument har vi undersökt den befintliga samlingen av probabilistiska integratörer för ODE, och föreslagit en ny konstruktion-den första att baseras på implicita metoder-ge en noggrann beskrivning av sina teoretiska egenskaper. Vi har gett preliminära experimentella resultat som visar effekten på parametern slutsatser om användningen av olika första ordningens metoder, både befintliga och nya, i utvärderingen av framåtmodellen. Överordnade flerstegsmetoder är tillåtna genom vår konstruktion. Vår diskussion om integratorkalibrering kräver inte någon lösning på detta subtila och taggiga problem, utan föreslår flera vägar för framtida forskning. Vi har diskuterat en fråga om förhållandet mellan skalningsparametern α och andra metodegenskaper. Insikt i denna fråga kan vara nyckeln till att göra dessa typer av randomiserade metoder mer praktiska, eftersom gemensamma trick för kalibrering kan uppstå som sedan är tillämpliga på olika problem. En intressant undersökningsriktning, som undersöks separat, handlar om huruvida uppskattningar av globala fel från andra källor, t.ex. Gränsöverskridande felmodellering, skattning av tillståndsnummer, skulle kunna användas för att kalibrera dessa metoder. Kommissionen ska anta delegerade akter i enlighet med artikel 10 i syfte att komplettera denna förordning genom att anta delegerade akter i enlighet med artikel 290 i fördraget om Europeiska unionens funktionssätt för att se till att bestämmelserna i denna förordning efterlevs.
Samtidigt som beräkningskostnaden för klassiska lösare multipliceras med mängden prover, kan dessa metoder fånga godtyckliga (icke-gaussiska) distributioner över lösningarna och kan minska övertro på omvänt problem för ODEs-som demonstreras i Conrad et al (2017, avsnitt 3.2). , Abdulle och Garegnani (2018, avsnitt 7) och REF.
Implicit Probabilistic Integrators for ODEs
{'venue': 'Advances in Neural Information Processing Systems 31 (2018) pp. 7244-7253', 'journal': None, 'mag_field_of_study': ['Computer Science', 'Mathematics']}